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Abstract 


Measurements from the space shuttle flights have revealed that a large spacecraft in a low 
earth orbit is accompanied by an extensive gas cloud which is primarily made up of water. The 
charge exchange between the water molecule and the ionospheric 0 + ions produces a water ion 
beam traversing downstream of the spacecraft. In this report we present results from a study on 
the generation of plasma waves by the interaction of the water ion beams with the ionospheric 
plasma. Since velocity distribution function is key to the understanding of the wave generation 
process, we have performed a test particle simulation to determine the nature of H 2 0 + ions 
velocity distribution function. The simulations show that at the time scales shorter than the ion- 
cyclotron period t c , the distribution function can be described by a beam. On the other hand, 
when the time scales are larger than r c , a ring distribution forms. A brief description of the 
linear instabilities driven by an ion beam streaming across a magnetic field in a plasma is 
presented. We have identified two types of instabilities occurring in low and high frequency 
bands; the low-frequency instability occurs over the frequency band from zero to about the lower- 
hybrid frequency for a sufficiently low beam density. As the beam density increases, the linear 
instability occurs at decreasing frequencies below the lower-hybrid frequency. The high frequency 
instability occurs near the electron cyclotron frequency and its harmonics. In the low earth orbit 
the low-frequency instability is likely to occur upto frequency / < f th = 5 kHz. The electron 
cyclotron frequency in the ionosphere is about 1 MHz ; thus an enhancement in the electric field 
noise upto a few MHz by the high frequency instability is expected. The wavelengths of the low 
frequency waves are found to range from a few tens of centimeters at relatively high plasma 
density (~ 1C $cm~ 2 ) to a few meters at relatively low densities (~ 10 4 c/n' J ) The wavelength of 
the high frequency waves is just a few centimeters. In order to study the nonlinear evolution of 
the waves, we have performed 2. 5 -dimensional particle-in-cell (PIC) simulations which show 
that transient amplitudes of the waves can be upto a few Volts / m at the lower-hybrid frequency. 
The frequency spectrum of the waves from the simulation shows that the water ion beam 
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produces a broadband electric field noise extending in frequency upto a few times the ion plasma 
frequency. The simulations also reveal an enhanced electric field noise level near the electron 
cyclotron harmonics. The spectrum is characterized by a very low frequency peak below the 
lower-hybrid frequency, a relatively broad peak near it, and peaks near the first few electron 
cyclotron harmonics. The high frequency waves have not been measured from the shuttle, 
probably due to the long antenna used in the measurements. 

2.5-dimensional simulations are computationally very expensive. Therefore, so far we 
have performed only a limited set of simulations which have not allowed us to reach a definite 
conclusion about the electric field levels in the different frequency bands. In order to reach a 
complete quantitative understanding of the effects of the contaminant ions on the space station's 
electromagnetic environment, we suggest that (1) further simulation studies be performed with 
realistic parameters of the ionosphere and (2) the high frequency response of the antenna used in 
the wave measurements aboard the space shuttle be further scrutinized to determine if the high 
frequency waves with wavelengths of a few tens of centimeters or shorter are artificially damped, 
causing a sharp roll-off in the frequency spectrum above about 10 kHz. 



1. Introduction 

Space shuttle flights have amply demonstrated that the contaminant gas plays an important 
role in determining the electromagnetic environment of a large spacecraft in the low earth orbit. 
In the case of the shuttle itself, water is the major contaminant [Paterson and Frank, 1 989], When 
water molecules undergo a charge-exchange with the oxygen ions, the major ion species in the 
ionosphere, a beam of H 2 0 + is created. Such ion beams traverse perpendicular to the 
geomagnetic field lines and interact with the ionospheric plasma to generate plasma waves, which 
contribute to the electromagnetic environment of the spacecraft. 

The space station, as planned, is a large spacecraft and its electromagnetic environment is 
likely to be affected by the contaminants. In this report we present results from a study on the 
contributions of the contaminants to the electromagnetic environment of the space station, using 
appropriate linear theory for plasma instabilities and their numerical simulations. 

A detailed summary of the measurements on the contaminant molecules and ions, and 
waves during some space shuttle flights can be found in a technical report by W. S. Kurth [ 1991], 
made to NASA Lewis Research Center. This report is a collection of papers in which 
measurements using the plasma diagnostic package (PDP) of the University of Iowa are reported 
and analyzed. Linear theories for the generation of plasma waves have been studied by Cairns and 
Gumett [1990, 1991], However, a linear instability theory cannot predict the saturated amplitude 
levels of the waves. In order to predict such amplitude levels, nonlinear theory and simulation 
studies are needed. The nonlinear processes also modify the frequency spectrum. 

Recent studies show that observed waves are best described as wave modes driven by a 
water ion beam traveling with the orbital velocity of 8 km/s , perpendicular to the earth's magnetic 
field, in the ionospheric O* plasma [Cairns and Gumett, 1991], The excited mode is lower- 
hybrid waves, which propagate nearly perpendicular to the magnetic field, and have a wavelength 
approximately given by X = 2 ni k ~ 2 / 0) o , where co 0 is the oxygen ion plasma 

frequency, and is the orbital velocity of the spacecraft. In low earth orbit the lower-hybrid 
frequency is about f fh =SkHz , and for the ionospheric plasma density in the range 
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10 10 - I0 12 m 3 , the expected wavelength may range from about 10 cm to a few meters. The low- 
frequency instability for / ^ f lh , driven by an ion beam traversing perpendicularly to the 
magnetic field, has been extensively studied in a laboratory plasma and a thorough discussion can 
be found in Seiler [1977], 

In this study, we first perform a test particle study to determine the nature of the water ion 
velocity distribution function. We show that only on a time scale of the order of the ion-cyclotron 
period ( r e ) the distribution function becomes a ring if the water ions are continuously produced. 
In the spacecraft frame of reference, in which contaminants are released into the ionosphere, the 
center velocity of the ring is equal and opposite to that of the spacecraft velocity, the minimum 
velocity is zero and the maximum velocity is double that of the spacecraft velocity. The ring is 
created by the cycloidal motion of the contaminant ions in the crossed magnetic field and the 
convection electric field. At shorter time scales t « r c , the ions are well described by a beam 

traversing perpendicular to the magnetic field in the ionospheric reference frame. We perform a 
linear instability analysis for the beams showing that growing waves are possible in two separate 
frequency bands, a low-frequency band ranging from zero frequency to the lower-hybrid 
frequency in the plasma. The high-frequency band covers the electron cyclotron frequency and 
its harmonics, which fall in the megahertz range. However, it turns out that for the ion beam 
velocity of 8 knvs, the growth rate for the excitation of the high frequency waves is about an 
order of magnitude smaller than that for the relatively low frequency waves near the lower-hybrid 
frequency. Therefore the high frequency waves may not grow to appreciable amplitudes. 

We also present results from numerical simulations on the nonlinear evolution of the 
instabilities. Using a 2. 5 -dimensional particle-in-cell (PIC) code, we have simulated the 
generation of waves by an ion beam in a background plasma. Simulations reveal the excitation of 
both the low and high frequency instabilities. However, due to the artificial electron to ion mass 
ratio, the ion beam velocity used in the simulations becomes unrealistically high relative to the 
electron thermal velocity. Therefore, the high frequency waves are excessively enhanced 
compared to that expected for the orbital velocity of 8 knvs in the ionosphere. However, the low 
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frequency instability is correctly simulated and simulations predict a broadband excitation 
extending well above the lower-hybrid frequency; considerably enhanced electric field noise is 
found to occur at frequencies several times the ion plasma frequency. This suggests wave 
generation in the frequency range of a few hundred kilohertz in the ionosphere. In addition, the 
high frequency instability causes wave generation in the megahertz range. PIC simulations with 
realistic parameters of the ionosphere are expensive to run. Therefore, we have performed a 
limited set of simulations which correctly describe the linear and nonlinear evolution of the waves, 
but at this stage the results cannot be directly applied to the ionosphere. In order to reach a 
definite quantitative understanding of the waves, we have recommended further simulation studies 
using the code developed here, and further scrutiny of the wave measurements aboard the space 
shuttle, especially for the high frequency waves with short wavelengths. 



2. Contaminant Ion Distribution Function 

Before a study on wave generation by contaminant ions is made, it is essential to know the 
nature of the ion velocity distribution function. The basic problem is as follows. As the neutral 
contaminant atoms and/or molecules are released in the ionosphere, they move with the orbital 
velocity of the spacecraft with respect to the ionospheric plasma. After they undergo a charge 
exchange reaction with the ionospheric 0 + ions, the contaminants become ions and they are 
affected by the ambient magnetic field and the electric fields around the spacecraft. Since in the 
low earth orbit, a spacecraft charges to a negative potential less than 1 volt, the electric field 
effects associated with the charging are expected to be insignificant, and the earth's magnetic field 
determines the nature of the velocity distribution function. We performed calculations on the 
trajectories of a large number of H-X)* ions when they are suddenly released in a magnetic field 
of B a = 0. 3 G . The initial velocity distribution function of the ion cluster was assumed to be a 
Maxwellian with a temperature of 300°K. In the rest frame of the spacecraft, the ions also see the 
convection electric field given by E = V_q x B^ , where is the spacecraft (orbital) velocity. 
Assuming B along the Z axis and along the -X axis, E is along the Y axis. The electric 
field causes an E x B^ drift in the X-direction, which is opposite to the velocity of the 
spacecraft. As the ions drift, their velocity distribution function changes in the V x -V y plane as 
shown in Figure 1, which gives the location of the ion cluster in the V x -V y plane at several 

selected times after their release. At the initial time / = 0, the average velocity of the ion cluster 
is zero. Due to the cyclotron motion, the trajectory of ions in the V x -V y plane is a circle given 


by 


F x = -f-(l + cos 0,1) 

B o 

E 

V v = — sin t 
y B n 


( 1 ) 

( 2 ) 


and the radius of the circle is U = E y I B a . In the above equations O w is the water ion cyclotron 
frequency. In our calculations we measured time in units of where fi 0 is the 0 + 

cyclotron frequency. The water ion cyclotron period is r cw — 2 7t I E2y/ or X2 0 
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= 2 nQ 0 / f2 v - 7.07. The normalized time appearing in Figure 1 is defined by t = ( Q 0 . For 
Figure 1, we assumed U = 1 5km i s. The trajectory of the ions in Figure 1 agrees well with that 
given by equations (1) and (2). We note that if the ions were free to execute their cycloidal 
motion, their velocity range wouid be from zero to 2 U in the spacecraft frame of reference. In 
the ionospheric reference frame the velocity ranges from -U to U over a cyclotron period, 
with an average velocity of zero. Near the release time, the relative velocity of the ion cluster in 
the ionospheric reference frame is -U , which is relevant for the wave generation process, which 
is much faster than the ion gyration. 

When the contaminant ions are created by charge exchange, they are created continuously. 
This continuous generation of ions modifies the distribution function significantly as shown in 
Figure 2. The ions progressively occupy a larger and larger portion of the circle shown in Figure 

1 and described by equations (1) and (2), and in about one cyclotron period, the circle is 
completely populated giving a ring distribution in the V x — V y plane. The thickness of the ring is 

determined by the initial thermal property of the ions. 

The above results are based on test particle calculations, in which waves generated by the 
streaming ions are not included. In a later section we describe results from a self-consistent 
simulation showing how the waves begin to mix the ions in velocity space at a time shorter than a 
cyclotron period. 
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3. Linear Instability Analysis 

We have identified two instabilities driven by an ion beam streaming perpendicular to the 
ambient magnetic field in a plasma. A low-frequency instability occurring for frequencies 
/ < f(h-> where f (h is the lower-hybrid frequency, and a high-frequency instability for 
/ = n f ce , where f ce is the electron-cyclotron frequency and n = 1,2,3,- •*. We briefly describe 
the linear aspects of the instabilities here. 

Low-Frequency Instability: Since the streaming ions excite waves at a time scale shorter than 
the ion-cyclotron period, it is useful to study linear wave excitation by an ion beam streaming 
perpendicular to the magnetic field. Assuming ions are cold and unmagnetized and electrons are 
highly magnetized, the dispersion relation for the angular frequencies a> «fi e = 2xf ce is given 
by 

1 + (l-a)<ng n) 

k 2 ( a-k x U ) 2 at * 2 at 

where co e and a> 0 are the electron plasma and O* ion plasma frequencies with the ambient 
plasma density N 0 , respectively; (o { is the contaminant ion plasma frequency with the density 
,V 0 ; a is the contaminant ion relative beam density, i.e., a = n w /N 0 with n, as the absolute 
density of the contaminant ions, k is the wave number with components k, and k x parallel 
and perpendicular to the ambient magnetic field, respectively; and U is the contaminant ion beam 
velocity along the x-direction. When the contaminant ions are H 2 0 + , cOj is denoted by 
o) w , and n, by n w . When rearranged equation (3) becomes 

Pat - 2fc x Pot + ( Pk\ - am Q - Q) + 2 k x Qa> -k x Q = 0 ( 4 ) 

where, 

P = 1 + y 2 cost 0, 

Q-\- a+Mcx sin 2 6 , 

y-(o e IQ e , k x = k x U / co 0 , oi-a}i co 0 , 

m 0 = m o l m i , m 0 g=m i /m e 
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where m g , rttj , and m e are the mass of an oxygen ion. a contaminant ion and an electron, 
respectively; and 0 is the angle between the wave vector k and a direction perpendicular to the 
magnetic field, i.e., k x = k cos# and k, - k sin 6. 

For small a and k z = 0, equation (3) can be solved analytically for the maximum 
growth rate and the corresponding wave frequency a> r and wave number k x . We seek solution 
for complex o) = oj r +iy and obtain 

k x = (1- a) 1/2 co th / U, co r = (1 - a) 1/2 co th (1 - S r ) , y= (1- a) 1/2 co (h 4 (5) 

where co (h is the lower-hybrid frequency, 8 r = 2~* /3 [ a / ( 1 - a)] l/3 and S t = J3 S r . These 

solutions are written by analogy to the dispersion relation for the two-stream instability 
[Hasegawa, 1975], 

We also solved equation (4) numerically and searched for complex roots a> = cb r +iy. 
The instability occurs when y > 0 . We have done this when the contaminant ions are water ions 
// 2 0 + . Figure 3a shows the variation of a> r = Re(<y) as a function of wave number k x for 
several values of the relative beam density a. The corresponding variation of the growth rate y 
is shown in Figure 3b. The plots in Figure 3a and 3b are for 6?=0, i.e., for propagation 
perpendicular to the ambient magnetic field. For a = 0.1, equation (5) yields results in 
reasonable agreement with a numerical results in Figures 3a and 3b. Note that for Figures 3a and 
3b we assumed Q e / co e = 1 / 3 . 

We notice from Figures 3a and 3b that growing waves occur for k x < 0.6, i.e., 
k x < 0. 6 coJU. The growing waves are the slow beam mode as indicated by the fact that all 
curves in Figure 3a fall below the slant line given by co = k x U or a> = k x . As the beam density 
increases, the frequency of growing waves is seen to decrease, i.e., the wave frequency roughly 
scales with the plasma frequency of the stationary ion population. The frequencies of maximum 
growth rate are indicated by dots, one on each curve in Figure 3a. The lower-hybrid frequency 
given by co, h = ^Jfi e Q = a > 0 / 3 is indicated by an arrow on the ordinate of Figure 3a. We see 
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from this figure that, depending on the relative beam density the frequency of the maximum 
growth rate can occur at frequencies widely differing from the lower-hybrid frequency. Only for a 
sufficiently small beam density does the frequency of the maximum growth (co max ) approaches 
Q)f h . For example, when the beam density is 10% (a = 0.1) and the stationary ion population is 
90%, the maximum growth rate occurs at co r ~ • 29 co 0 ~ (o fh . As the beam becomes denser, 
< 0 ,,^ decreases. When a = 0.9, =.09o) 0 = 0.2S(Of h . Figure 3b shows that the maximum 

growth rate is obtained when a = 0.5; as the relative beam density departs from this value, 
whether it increases or decreases, the growth rate decreases. It is also worth mentioning that the 
growth rates are approximately the same when the relative beam density is a or 1 - a 

Figure 3 b shows that the growing wave number is limited to k max <, 0.6 or 

k x <, 0.6 o) 0 l U . For the maximum plasma density of \&cm~* in the ionosphere, the maximum 

value of <a 0 = 327 krod/s. For this maximum value of co 0 along with U = & km / s, 

k x <, 24 m~ l or the wavelength X > 2n!k = 25cm. If the plasma density is lower, the 

wavelength becomes longer. For example if n 0 = 10 4 cm -3 , X > 250cm. Since the ionospheric 
plasma density typically ranges from 10 4 to 10 6 cm -3 , the wavelengths of the linear waves can 
range from a few tens of centimeters to a few meters. 

Figures 3a and 3b are for 0- 0. The variations of the wave frequency and the growth 
rate with 0 are shown in Figures 4a and 4b, respectively. In these figures the curves are the plots 
of the frequency ( (b r ) and the growth rate (/) at wave numbers where the latter ( y) maximizes. 
Curves for #>.01 are not reliable because at such angles the Landau damping by electrons 
reduces the growth rates and we have ignored this damping in writing the dispersion relation in 
equation (3). It is worth noting that the range of frequency for the wave excitation is significantly 
enhanced when d> 0. For example, when a = 0.1, the growing waves occur upto co r = 0.6 co 0 
for the angles 6 ^ 0.01 rad compared to ty r =0.36<u o at d- 0. 

High Frequency Instability: Besides the low frequency waves below the lower-hybrid 

frequency, there Is a possibility of exciting high frequency waves at frequencies near / = nf ce , 
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i.e., the harmonics of the eiectron-cyciotron frequency [Forslund et aJ, 1970], The excited waves 


are the electron cyclotron waves described by the dispersion relation 


{k x X d )i = -\ + e--I 0 {X) + 2ar-f 


SI 2 ( k x V w } 2 Vo 


( 6 ) 


where A = (kr e ) 1 /2, r e is the electron Larmor radius given by r e = V e If2 e , /„(•) is the 
modified Bessel function of the first kind, X d is the plasma Debye length, and the terms 
containing Z' represent the contribution of ions to the dielectric function of the plasma. Z' is the 
first derivative of the plasma dispersion function [Fried and Conte, 1961], The fourth term on the 
right hand side of equation (6) is the contribution from H 2 0+ ions and the argument of Z' in 
this term contains the Doppler-shifted frequency co-k x U divided by the product of k x and the 
thermal velocity of the water ions V w . Likewise, the fifth term is the contribution of oxygen ions. 
In the limit of small drift velocity, the real frequencies co r and growth rates y of the waves are 
given by [Forslund et al, 1970] 


co r =nQ e 


r = 


nQ 1 


ZIMr) 


2^2 ^ [\+(kx d ) 2 -z;(4 r )f 


where 


4 r = (a> r -kU)/kV w , Zl m () = Im(Z'(-))> Z;() = Re(Z'(» 


(7) 

( 8 ) 


(9) 


The growth rate maximizes when Z.'_( ) maximizes. Since Z’ m (% r ) = - 2£ r rt n e 




which becomes maximum when % r = -\l*j2 and the maximum value of Z^t) = 1.5. The 


corresponding value of Z'(g r ) = 0.53. Combining these facts, we obtain 

y _ 1.5n 1 1 


( 10 ) 


Q e 2rt a kr e [l.S + ^/lrf) 2 } 2 

The corresponding wave number k can be calculated from (co- kU) / kV w = - 1/V2 or 
co = k(U + V W / >/2). Since U » V w , k = rtO e / U or kr e = nV e / U , and the growth grate is 
given by 
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1 


( 11 ) 



1.5 U 

2/r" 2 V e ( . ; Q e V} ' 

(1 ' 5+M 


Since V e »U ( U=%km/s , and 1' = 1 90£m / 5 for O.leV electron temperature), 



rc 4 


1.5 

2^ /2 


(^) 5 (^) 4 


a 


( 12 ) 


From the foregoing expression for the growth rate it is clear that the wave excitation near the 

electron cyclotron frequency and its harmonics is favored by a large plasma density 
so that o) e »Q . Assuming typical values of U - Skm / s and V e = 190 brt/s, 
Y I = 5.3 x lO" 8 /?- 4 ^) 4 In a relatively dense ionospheric plasma with a density of 10 6 cm -3 , 

co e = 2/rx 9 x \Qfiradis, the cyclotron frequency = 2;rx8.7 x 1C Prad/s and the 
corresponding growth rate r = 6 x 1 0^ Q e . 

We found that the growth rate at low frequencies {f < f Q ) is about 0.1 co 0 , which is 
about an order of magnitude or more larger than the corresponding growth rate at the frequencies 
near the electron cyclotron frequency and its harmonics. Therefore, the low frequency waves 
should dominate in the ionosphere. 
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4. Numerical Simulation 

Simulation Technique: In order to estimate the effects of noniineanties on the frequency 

spectrum and the wave amplitudes, we have performed numerical simulations using a particle-in- 
cell (PIC) code. We performed 2. 5 -dimensional simulations as illustrated in Figure 5. The plasma 
is two-dimensional in the x-y plane; the ambient magnetic field is primarily along the : axis, but it 
has a tiny component in the x-y plane as well. The plasma particles are like rods of finite size 

[e.g., see Birdsall and Langdon, 1985]. Although the simulation is two-dimensional in the 
configuration space, in velocity space we include all the three velocity components, i.e., V x , V y 

and V z , For this reason the simulation is called 2.5-dimensional. We solve the equation of 
motion of charged particles in given magnetic fields and electric fields determined by the solution 


of the Poisson equation given by 

(Pep <p 

+— = -pi £o 


(13) 


dx2 dyl 

where <p is the electric potential, s a is the permittivity of free space, and p is the charge 
density determined by the plasma particles. We solve equation (13) subject to the periodic 
boundary conditions in both x and y directions. The numerical solution of equation (13) is 
implemented by the Fast-Fourier-Transform (FFT) method, a two-dimensional Fourier transform 
of equation (13) gives 


~<KK,k y ) = 


1 pQCxJCy) 

«„<*!+*)) 


(14) 


where k x and k v are the wave numbers associated with the variation of the potential along the 

x and y directions, respectively. First we perform a two-dimensional transform of 
p(x,y) to obtain p(^ x ,^ v ) using a Fast-Fourier-Transform IMSL subroutine. The solution 

0(x,y) from <fi(k x ,k v ) is obtained by an inverse transform using the same subroutine. 

The equation of motion of charged particles is given by 

dV„ 


m a -jY = q a {E + V a <B^) 


(15) 
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where m a and q a are the mass and charge of a charged panicle of type a , V a is its velocity, 

E_ is the electric field and is the ambient magnetic field. In our model (Figure 5), 

E = a x E x + a y E y , B^ = B x0 a x + a y + B zo a z , and equation ( 1 5) can be written as 


m n 


dV, 


dt 


m a 

m a 


dV, 


va 


dt 

dV ta 

dt 


= q a (E x+ V v B zo -KB^) 

= q a ( E y - V x B zo + E z B xo) 

^qaiKByo-VyBxo) 


(16) 

(17) 

(18) 


Equations (16) and (17) are advanced in time using the leap-frog technique [Morse, 1970], 


Equation (18) is advanced by realizing that 

-77 ~ <la * B vo + <la y B xo ) = 0 

u t 


(19) 


or 

m^za -q a xV yo+<la y B xo = constant (20) 

With the known initial velocity and position of a charged particle, its velocity component at later 
times is determined in terms of its x(t) and y(t) coordinates which are available by solving 

(16) and (17). 

We use the following normalizations and definitions in our subsequent discussions: 
position r=(r/A rfo ), time t=tco e , velocity V = F/I7, density N=N!N 0 , potential 
^ = <f > / <f > n , and electric field £ = E/E Q , where X do , co e ,V e and N a are the Debye length, 
electron plasma frequency, electron thermal velocity and density in the initial unperturbed plasma, 
respectively; <f>„= k B T a ! e , E a = <p n ! X do , and T 0 is the initial electron temperature in the 

plasma. 

Figure 5 shows that the simulation system is gridded with grid spacings Ax and Ay in 
the x and y directions, respectively. The system size L x x Ly and the grid spacings 

Ax and Ay were varied during the course of the simulation. However, a typical system size is 
L x x Ly - 256 x 128 with Ax = 2 and Ay = 2 . 
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Initially at time t = 0, total numbers of electrons and ions were loaded uniformly in the 
x-y plane. The total number ( N p ) of electrons and ions in the system are equal and typically 

-Vp = 163840, which gives 5 particles per Debye square. In the simulations, each particle 

represents a large number of real electrons or ions. We used two types of ions corresponding to 
the ambient O* ions in the ionosphere and the contaminant ions such as the water ions, which 

are produced by charge exchange between 0 + and H 2 0 . In the rest frame of the spacecraft, at 

the time of the charge exchange H 2 0 + ions are stationary while 0 + and electrons drift relative 

to the spacecraft. This relative motion generates waves, which we are trying to simulate here. 
We performed our simulation in the ionosphere reference frame, i.e., we assume that O* and 

electrons are stationary and H 2 0 + are moving. 

The initial velocity distribution functions of the electrons and ions are given by 


F ' = 


Fw - T T — 7 7, m w W,-U?+V} + V})lk B T c \ 

(2xk B T 0 lm w y a 1 y 




where T a is the initial temperature of electrons, oxygen ions and water ions; m e , m 0 , and 


are the corresponding masses of the charged particles, U is the drift velocity in the rest frame of 
the ionosphere, with respect to which the H 2 0 + ions drift. In the above expressions, k B is the 


Boltzmann constant, N a is the unperturbed density and the ion densities n 0 and n w are such 
that N 0 = n 0 +n w . 


Two-dimensional simulations with real masses of electrons and ions become prohibitively 
expensive. Therefore, we have performed simulations using artificial mass ratio; we have 
performed simulations with m 0 / m e = 5 1 2 and 2048 while keeping m 0 / m w = 1 6 / 1 8 . The last 

ratio correctly reflects the ratio of oxygen ion mass to the water ion mass. The former ratios are 
chosen to be sufficiently large to adequately separate the electron and ion time scales. As we 
discuss later, the larger the mass ratio m Q / m e the more realistic the simulation results become. 
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Typically in the ionosphere 0+ thermal velocity is about V 0 = \ km/ s, electron thermal 

velocity V e = 190 /km / s, and the orbital velocity about U =Sknt / s. For the above reduced ion 
mass in the simulation, the thermal velocity of 0 + ions V 0 = ' 2 V te , which is about an order 

of magnitude (or more) larger than the real value in space which is V Q = 1/ 190 V e = 0.0053F e . 
Since for the low-frequency waves the ion dynamics is important in the present study, we keep the 
ion thermal and drift velocities in the same proportion as in the ionosphere, but then the drift 
velocity becomes a sizable fraction of the electron thermal velocity. For example, we present 
results from simulations with m Q l m e = 512 and U = 9V o ^0AV e , and also for much heavier 
ions with m a /m e = 2048 and U = 0.2V e . These choices of the drift velocity for the respective 
mass ratio correctly predict the low frequency behavior, but they also reveal enhanced electric 
field noise near the electron cyclotron frequency and its harmonics. 


Numerical Results: 

Before we begin to present numerical results, it is worthwhile to know some of the 
characteristic plasma frequencies in low earth orbit. Theses frequencies are summarized here as 
follows: 


electron cyclotron frequency: 


fee 


1 eB o 


m„ 


0.87 MHz 


0 + cyclotron frequency: 


fa 




~~ fee = 30 //r 


lower-hybrid frequency: 
electron plasma frequency: 


f th =4UZ = ^Hz 


n n e A All 


/ e = U(^) S1-9M: 


m e s 0 


ion plasma frequency: f Q - 5.2-52 kHz 

The cyclotron and the lower-hybrid frequencies are for a magnetic field B 0 - 0.3G. The 
values of f e and f Q are based on plasma densities in the range N a = 10 4 - 10 6 cm 3 
Comparing f ee and f e , we notice that f ce / f e < 1 ■ In the simulation we typically chose 
fee I fe~\ 0.3, which represent the situation in the wake region and outside it, respectively. 
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In our calculations lower-hybrid frequency is an important parameter; the value quoted 
above is approximate and its exact value is given by 

fa, = JZZ, ( , (24) 

When f ct f ci « / o 2 , f lh = ajfgJa . On the other hand when the plasma density becomes 
sufficiently low making fcefci » fo » flh = fo‘> ' n ionosphere this should happen when 
plasma density N a < 10 4 cm -3 and it is a distinct possibility inside the wake of a spacecraft. 

With the artificial mass ratio m 0 lm e - 512, the ratios of some of the frequencies given 
above are modified as follows: f e l f Q - 22.6 and not 171.4 like in the ionosphere, and similarly 
fee 1 f co = 512 and not 2. 9 x 1 0 4 . Since ions play a crucial role in the low-frequency instability 
described here, the lowering of these ratios allows us to see the complete evolution of the 
instability in a comparatively short time. The lower-hybrid frequency in our simulation is given by 



For a typical ionospheric density of N a = lC^cm -3 , B 0 = 0.3 G and m 0 !m e - 512, f ce I f e = 0.3 

and =0.013 f e . When the mass ratio m 0 ! m e is increased to 2048, f (h = 0.067 f e . When 

the plasma density is reduced, like in the wake region of a spacecraft, the ratio f (h / f e is 
increased; for example if N a = 1 0 4 c/n -3 , f ce i f c = 1 and fff, ! f e = Q 03 for m 0 im e - 512. 

Simuiatiou with m 0 lm e - 512, U/V e = 0.4 and f ce ! f e - 1: 

This simulation was performed to check some of the linear features of the waves as 
described in section 3. We performed the simulation for f ce l f e = 1 , which is the case in the very 
low density region of the wake. As noted above, this gives f th = 0.03 f e . We performed 
simulations with I N o =0.\, T 0 -T w -T e -T 0 and U / V e = 0.4. One may wonder why such 

large drifts in terms of electron thermal velocity are chosen. As explained earlier, tor the artificial 
mass ratio m 0 im e - 512, the ratio of 0 + thermal velocity to that of the electrons is 
V 0 !V e =0.045 and not 0.006 as for the real mass ratio of m Q I m e = 29376. In order to 
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properly treat the low frequency waves, for which ion dynamics is crucial, we keep the ion Mach 
number, i.e., M — U !V 0 as close as possible to that in low earth orbit. The choice ot 
U / V e = 0.4 gives AT = 9 which is close to that given by a drift velocity of 8 km/s. 

Figure 6 shows the temporal evolution of the wave; in Figure 6a we have plotted the 
normalized electric potential ~<p at the center of the simulation region, i.e., at the point (128,64). 
Figures 6b and 6c show the electric fields £ x and £ y at the same point, respectively. These 

figures clearly show that the waves grow and then they decay. It is also apparent that the waves 
have low and high frequency components. From Figure 6a, we find that the dominant low 
frequency component has a time period r f =260o)~ l , which gives a frequency 

/ = t; 1 = (2;r/260)/ e = 0.024/ e , which is slightly lower than the lower-hybrid frequency 

mentioned above. We saw in section 2 that for sufficiently low beam density, the strongest 

growing waves have frequencies close to the lower-hybrid frequency (see Figures 3 a and 3 b). 

The maximum amplitude of the electric field is about 0.25kT o / eX d . For kT 0 = Q.2eV 
and X d -ZAcm for a plasma density N o = \0*cm~ 3 , the maximum electric field 

E = (E 2 + E 2 ) V2 « >/2 x 0.25x0. 2/. 034 = 2V /m. Thus in response to the injection of water and 

x y 

its subsequent charge exchange, transient fields of a few Volts/meter are expected, which is in 
contrast to the average rms fields of about a few multivolts/m measured from the shuttle [Gumett 
et al, 1988], 

From the electric potential and the fields, one can estimate the wave number and 
wavelength of the waves. Since E — kjj>, where k = k x a x + k v a y , we estimate that a typical 

value of k x for the low frequency wave is k x - £ x I~<j>~0 1, which yields a wavelength 
X x = = (2sr/0.0)X d = 89 X d . The linear solution for the waves in equation (5) yields 

k x = 0. 066 X~J and X x = 95X d . Thus the wave number from the simulation is seen to agree with 

the linear solution. 

The frequency spectrum of the electric fields shown in Figure 6 is plotted in Figure 7. 
There are two curves in this figure. The solid curve is obtained by performing Fourier transform 
of the data plotted in Figures 6b and 6c; the data analyzed is obtained by calculating 
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E = (E 2 ^ E ~) V2 . The broken-line curve is the noise in the simulation when the drift velocity 

U = 0. We notice that the spectrum is considerably enhanced over the noise level in two 
frequency bands; the low frequency band is f <0.1 f e , which in our simulation is about 
2 f Q , where f a is the ion plasma frequency. The high frequency band consists of the sharp 
peaks near the electron cyclotron harmonics. These enhancements are expected from the linear 
instability analysis presented in section 2. 

The low frequency band has its own structures. We see a strong enhancement near the 
lower-hybrid frequency which in our simulation is 0.03 f e . We also see a strong enhancement at 
very low frequencies f <0.01 f e . The peak near the lower-hybrid frequency is expected from the 
linear instability analysis presented in section 3; but we see that it is considerably broadened, and 
the enhancement extends to higher frequencies upto twice the ion-plasma frequency. The 
enhancement at very low frequency represents the temporal modulation of the wave amplitude 
shown in Figures 6a, 6b, and 6c; this modulation consists of the initial stage of the wave growth 
followed by its saturation and eventual decay. 

In the simulation described above we have used a relatively low artificial mass ratio of 
m 0 ! m e - 512, which reduces the ion cyclotron period. Thus, as the ions excite the wave they are 

significantly affected by their cyclotron motion. This is shown in Figure 8, which displays the 
evolution of the ions in the V x - V v velocity space as the instability develops (Figure 6). Water 

ions are injected at / = 0 with a drift velocity U = 0.4V e along the x direcuon. 

The ion cyclotron period in the simulation for water ions is 
t c = 2kI 0^ = 27^171^1 m e )l Q e . Since f2 t = a) pe in the simulation described here, 

r c = t c co pe = = 3619 . This cyclotron period is clearly visible from the rotation of 

the water ions around the oxygen ions in the V x - V v plane. We saw a similar rotation in Figure 

1. The cyclotron motion shown in Figure 8 modulates the waves contributing to the very low 
frequency peak in the frequency spectrum in Figure 7. In reality the normalized cyclotron period 
is considerable larger and the mixing of water and oxygen ions as seen in Figure 8 for t > 2000 
occurs at a much shorter time scale relative to the cyclotron period. Thus the modulation due to 
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the cyclotron motion may not produce as low frequency as the ion cyclotron frequency, but a low 
frequency component is expected. 

Simulation with m 0 l m e = 2048, U / V e = 0.2 and f ce / /, = 0.3: 

The simulation with the mass ratio of 2048 is more realistic, but we are still away from the 
real mass ratio by a factor of 14.3. Due to constraints on computer resources, so far we have not 
done the simulation with the real mass. Figure 9 shows the temporal evolution of the waves seen 
in the simulation. In this simulation, we used n o = n w =0.5N o . Like in Figure 6, Figure 9a 

shows the evolution of the potential at the center point of the simulation region, while Figures 9b 
and 9c show the evolution of E x and E v at the same point. Note that the wave is seen to 

grow to a maximum amplitude of j^j = 4 kT 0 / e and the corresponding electric fields E x and E y 

have amplitudes upto about 0A5kT o /eA d . For kT a = 0.2eV and A d = \cm, the maximum 
electric field seen from this simulation is E « (E 2 + E 2 ) l/2 = -Jl x0.15x.2/0.01 = 5V/m, which 

is somewhat larger than that seen from the previous simulation. The frequency spectrum of the 
waves shown in Figure 9 is plotted in Figure 10. Note that in this simulation f ce ! f e = 0.3 and 
the 0+ ion plasma frequency is f Q - 0. 022 /, . The broken-line curve in Figure 10 is the 

background noise of the simulation when U - 0. The solid-line curve gives the frequency 
spectrum of the electric field data for E - ( E ' 2 M E~) l/2 shown in Figure 9. The most dominant 

peak in the spectrum occurs below the lower-hybrid frequency. The enhancement in wave power 
over the background noise by an order of magnitude or more occurs upto a frequency / =0.1 f e 
which is about 5 f Q . Figure 10 also shows the high frequency waves near 
/~0 3/, , 0.6 /, and 0 9/,, which are the harmonics of the cyclotron frequency f ce =0.3/,. 
The evolution of the ion phase space in the V x - V v plane is shown in Figure 11. In the 

present simulation which we are discussing here, the ion cyclotron period ~r c - 42892 . Figure 1 1 
shows that the water ions with the initial drift velocity of U = 0.2F„ begin to gyrate as seen 
from the plots at / = 0, 2000 and 4000. However, the waves act on them and they mix with 
the 0 + ions. The mixing has begun at t = 4000 and by the end of the simulation at f = 9000 
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the two types of ions have undergone a considerable mixing in the l x V y velocity space. This 
feature of the present simulation with mass ratio of m Q i m e — 2048 and f ce l f e - 0.3 is quite 
different from that shown in Figure 8 from the simulation with m a i m e = 5 12 and f ce / /* = 1 
The mixing of the ions in the X - V x phase space is shown in Figure 12. At the initial time we 
note the O* ions have zero average drift velocity while the water ions have an average drift 
U = 0.2. Figure 12a shows that at t =1500, there is only a weak perturbation in the phase 
space. As the waves intensify (Figure 9), the two ion populations undergo a coupling through the 
electric fields of the wave. As seen from Figure 12b, vortices appear by the time J = 3500 and 
they fill up progressively mixing the 0 + and HXX ions. 

The ion distribution in the V x — V y plane shown in Figure Ilf is a pancake type as 

observed from the space shuttle [Paterson et al, 1989], However, the distribution shown in this 
figure is in the frame of the reference of the ionospheric plasma in which O* ions are stationary 
while the FIX)* ions drift with the shuttle speed. The waves modify the distribution functions of 
both oxygen and water ions in the plane perpendicular to the magnetic field. Figures 13a and 13b 
show the parallel velocity distribution function F(V X ) of O* and FIX)* , respectively Over the 
simulation time of I = 9000 , O* ions show hardly any change in F(V Z ), while HXX ions 
show a slight modification; these ions gain a slight drift along the z direction due to the V y x B x 
force acting on them. Figures 13c and 13d show F(V X ) and F(V y ) for the O* ions, 

respectively. The O* ions are seen to develop nearly a symmetric tail in V x which amounts to 
the ion heating. The O* ions are seen to develop an asymmetric tail along X y . This is 

connected with the cyclotron motion of HXX ions and the interaction between O* and HXJ* 
ions through the waves. Figures 13e and 13f show the F(V X ) and F(V y ) for the ions, 

respectively. It is seen that the FIX X ions loose their initial drift along X to that along V y due 

to the cyclotron motion. It should be also pointed out that partly the loss of the momentum, i.e., 
reduction in V x must be due to the wave generation. Besides the transfer of drift motion from 
V x to V v , it is also seen that the water ions are heated as evidenced by the increase in the 

velocity spread of the distribution shown in Figures 13e and 13f with increasing time. 
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The electron distribution functions corresponding to those of ions in Figure 1 3 are shown 
in Figure 14. Figures 14a, 14b and 14c show F(V X ), F(V v ) and F(V Z ), respectively. These 

figures shown only slight heating in V x and V y . On the other hand, the heating in V z is 

appreciable. The electron and ion heatings seen here are similar to those reported by McBride 
[1972], who studied the generation of the lower-hybrid waves in a plasma where the entire ion 
population drifted across the ambient magnetic field. 



5. Conclusions and Discussions 

We have performed a systematic study of the wave generation by contaminant ions near a 
large spacecraft in low earth orbit. Our study includes (1) test particle simulation on the 
contaminant ion velocity distribution function, (2) linear wave analysis, and (3) self-consistent 
numerical simulation of the wave generation by the contaminant ion beam traversing 
perpendicular to the electric field. 

The test particle simulation revealed that for short time scales at which the waves are 
generated the contaminant ions can very well be represented by an ion beam rather than a ring, 
which forms at the time scale of the ion cyclotron period. The self-consistent 2.5-D numerical 
simulations with large ion to electron mass ratio corroborate the above statement. 

The linear instability analysis revealed that the waves occur in two frequency bands: ( 1 ) a 
low frequency band for f <f (lt , and (2) a high frequency band with / = n f ce , where f ce is 
the electron cyclotron frequency and n - 1,2,3, ••-. Comparing the growth rates of the waves in 
the two frequency bands for a contaminant ion beam velocity of U - &km / s, it turns out that 
the waves in the high frequency band may not be appreciably amplified because U «V te , the 
electron thermal velocity. In view of our simulation with a relatively large mass ratio, for which 
the high frequency waves remain significantly amplified (Figure 10), the amplification needs to be 
firmly established by further simulations as we discuss later. 

The self-consistent 2.5-D simulations do show that the lower-hybrid waves are excited 
along with the high frequency waves. The low frequency waves can have wave amplitudes upto 
a few Volts/m depending on the relative densities of the contaminant and oxygen ions. 
Simulations show that the wavelength of the low frequency waves is in agreement with that found 
from the linear instability analysis. 

In contrast to the linear-instability analysis, the 2.5-D simulations show low-frequency 
wave excitation in a broad frequency band, not limited to / < f (h , but extending to frequencies 
several times the Q + ion plasma frequency. The extension of the low frequency waves above the 
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lower-hybrid frequency is not predicted by the linear instability analysis, it is a consequence of 
nonlinear plasma processes. 

The low frequency waves at frequencies f < 30 kHz are not relevant to the space station s 

electromagnetic environment. In low earth orbit the lower-hybrid frequency is about 5 kHz and 
the oxygen plasma frequency f Q can range from 5 to 50 kHz as the plasma density ranges 
from 10 4 to lO^m -3 . If the nonlineariy generated low frequency waves extend in frequency 
upto several times f Q , it is distinctly possible that waves in the frequency range upto a few 
hundred kHz are generated. The observations from the space shuttle show waves extending to 
50 kHz. However, the measured frequency spectrum (see Figure 1 5) shows that the wave power 
sharpiy fails off at frequencies above the lower-hybrid frequency [Cairns et al, 1991]. It is worth 
pointing out that with increasing wave frequency the wavelength get shorter and shorter, making 
wave measurements with a relatively long antenna rather inaccurate; the antenna response may 
a tt pnn 3 t<* the measured waves. The IOWA PDP instrumentation carried a 3.9-m long antenna 
[Gumett et al, 1988], which could not correctly measure waves having wavelengths shorter than a 
few tens of centimeters. 

The comments made above also hold good for the high frequency waves near f - n f ce . 
The wavelengths at such frequencies are shorter than the electron gyro-radius which is about 3 
cm in low earth orbit. Waves with such short wavelengths are difficult to measure with an 
antenna of a few meters in length. Comparing Figures 10 and 15, we note that the simulated and 
measured spectra have some similarities, at least for frequencies / <30 kHz. The difference at 
higher frequencies is due to simulational inadequacies or due to the error in the measurements. At 
this juncture we are not sure about the real cause of the difference. 


Recommendations: 

Under the present grant, we have developed the analytic capability to predict the entire 
frequency spectrum of the waves generated. The main tool we developed is the 2.5-dimensional 
PIC code. So far we have done simulations with only 5 particles per Debye cell, which is too 
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small and it gives enhanced plasma fluctuations even without the wave generation. The 
background noise in Figures 7 and 10 are the consequence of this limited number of particles in 
the simulation. So far our simulation is also limited in terms of ion to electron mass ratio; we 
have used m 0 l m e = 2048 , which is smaller than the real mass ratio of m 0 / m e by a factor of 
about 14.3. The use of artificial mass ratio is required for the purpose of the computational 
economy. The run with the mass ratio 2048 used 18 hours of CPU on a VAX 7000 machine. 
The use of the artificial mass ratio also necessitates that the drift velocity of the water ions be 
relatively large in terms of electron thermal velocity. This enhances the wave generation at the 
electron cyclotron frequency and its harmonics. 

We recommend that a few simulations be performed with the real parameters of the 
problem, and with a sufficiently large number of particles that the background noise is sufficiently 
suppressed. Such simulations along with further scrutiny of the measured wave data in terms of 
antenna response for correctly measuring the short wavelength (high frequency) waves are 
essential for a full understanding of the effects of the contaminant ions on the electromagnetic 
wave environment of the space station. 

The low frequency waves with frequencies upto a few tens of kilohertz do not appear to 
be of much concern for the space station. But the high frequency waves in the MHz range near 
the electron cyclotron frequency and its harmonics may be of concern. These waves are 
extremely short wavelength waves. We feel that the PDP measurements at such high frequencies 
may not be correct due to the antenna response. Further simulations, as suggested above, may 
prove to be useful. 
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Figure 2. Ions' trajectory in the V x - V v plane when they are continuously injected, (a) 
/: 2i^ ! ,(b) t - and (c) / = 8/T^ 1 . Note that with continuing injection of ions the beam 

arc at early times evolves in a ring distribution. 
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Figure 3. Results from linear instability analysis, (a) Wave frequency co r R? (oj ) as a function 

of wave number k x . Note that 0) o is in units of the oxygen ion plasma frequency determined by 

the density in the ambient plasma without the water ions, (b) Growth rate is in units of co 0 ■ ^e 
several curves m (a) and (b) are for different values of the water ions relative density given by a 
in the legend. Other parameters for this figure are 9=0 and co e /fl e = 3 . 







Figure 5. Geometry of the 2 1 /2-dimensional simulation. Simulation is performed in the x-y 
plane. The dominant component of the magnetic field is in the z direction. A tiny component o 
the magnetic field is in the x-y plane; in the simulations described here this component is ong 
the x direction. 
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Figure 6. Temporal evolution of (a) normalized potential 0, (b) normalized electncjield 
E x and (c) E v . The horizontal axis is normalized time t =ta) e . For this run f2 e -co e , 

n=0.\N o , m,im e = 512, and U/V e = 0.4. 







Figure 7. The frequency spectrum of the electrical field E = (E 2 + E 2 ) l/2 is shown by the solid 
line. The electric field data for E x and E v are the same as plotted in Figures 6b and 6c. The 
broken line shows the frequency spectrum with zero drift for the water ions. 




Figure 8. The gyration of water ions in the V x - V v plane is shown for co e -O e . Eventually the 
water and oxygen ions mix the velocity space as seen for t > 2000 . 








Figure 10. The same as Figure 7, but for Q e - 

U = 0.2V e . The frequencies shown on t 






Figure 11. The same as Figure 8, but for O e = 0.2co e , m 0 ! m e - 2048, and U -0.2Y 
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Figure 13. Ion velocity distribution functions from the smulation wit ' A 

A AT. -0.3. m o^ m e ~ 2048, and f/=0.2F <•> A' AmA £ . fhTsodi 

F(K,) for O*. (d) F(V v ) for O'. (e) F(V X ) for Hfi , and (0W,) - 

curves are for f = 0, dasher curves for ,' = 6000. and .he dashed curves for < =9000. 
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Figure 14. Electron distribution functions from the run with Q e - 0.3 co e , n w /N 0 - 0.5, 
m 0 l m e = 2048, and U = 0.2V e . (a) F(V X ), (b) F(V y ) and (c) F(V Z ). The solid curves are 

for t = 0 and the dash-dot curves for t = 9000 . 
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Figure 15. Typical frequency spectra of the electrostatic broadband noise measured from the 
space shuttle fKurth, 1991; Cairns and Gumett, 1992]. 
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c* 

c* 

c* 


2.5 Dimensional Particle-In-Cell 
by Contaminant Ions Near a large 


Code For Studing Wave Generation* 

spacecraft. * 

* 


C* Theta = 0.05 (Angle of the Magnetic field to the X-axis) 

C* Bm = 0.3 (Normalized Magnetic field) 

C* Bxm = Bm*Sin( theta ) 

C* Bym = 0 

C* Bzm = Bm*Cos ( theta ) 

C* dx = 2 (Delta x) 

C* dy = 2 (Delta y) 

C* Time = Normalized time 


c* 

dt = 

0.2 (Delta Time) 


* 

c* 

Mass 

ratio of Oxygen Ions 

= 2048 

* 

c* 

Mass 

ratio of Water Ions = 

2304 

* 

c* 

Only 

Water Ions have drift 

velocity of 0.2 

* 


C* Number of grid in x-direction = 128 * 

C* Number of grid in y-direction =64 * 


REAL ION 

COMMON/PAR/ELE( 163840 , 5), ION (163840, 5) 

COMMON/ARRAY 1/PH I (129,65) 

COMMON/EFI ELD/EX ( 129, 65), EY (129, 65) 

COMMON/ARRAY/RHO ( 129,65) ,RHOE( 129,65) ,RHOI (129,65) 
COMMON/DEL/DVXI , VXIW ( 301 ) , FIW( 301 ) ,GIW( 301 ) ,VDRI ,DVXE,VDRE 

> ,VXIO( 301 ) , FIO( 301 ) ,GIO( 301 ) 

> , VXEO( 241 ) , FEO( 241 ) ,GEO( 241 ) 

COMMON/VEX/XX ( 1 2 9 ) , YY ( 6 5 ) 

COMMON TIME , NI J 

COMMON NX , NY , NXl , NYl , P I , MAIN 

COMMON RATO , RATW , DELT , VELVAR , BXM , BYM , BZM , EYO 
COMMON I SEED , DY , DX ,NIED,NIEDl ,NSUMO 

DELT=0 . 2 
TSTART = 0 . 

TMAX = 9000. 

I TST = INT( TSTART) 

NIT = ( TMAX- TSTART ) /DELT + 1 

NX=129 

NY=6 5 

NXl=NX-l 


NYl =NY- 1 


NIJ = 20 
I SEED=1 0 0 0 00 01 


c 


CALL SET 
CALL DISTELE 
CALL DISTIONW 
CALL DISTIONO 
CALL BEGIN 
CALL INITIAL 

TIME=-0 . 2 + TSTART 
LWRTE = 1 
LAC = -1 

OPEN ( UNI T = 9 2 , F I LE= ' EXEY-MASSO10Wl0Fixd0.DAT' , STATUS= ' NEW' ) 
OPEN( UNIT = 9 3 , FI LE= ' PH I -MAS SOI 0W1 0 FIXD0 . DAT ' , STATUS= ' NEW' ) 
OPEN ( UNI T = 9 4 , FI LE= ' RH O-MAS SOI 0W10FIXD0 . DAT ' , STATUS= ' NEW' ) 

DO 1 MAIN=0 , NIT 
TIME=TIME+DELT 
LAC = LAC + 1 
CALL DENS 
CALL POTENT 
CALL FIELD 



c 

c 

c 

c 


41 

42 
C 
C 


c 

1 

c 


c*** 

c*** 


c 


c 


19 

55 

C 


17 

56 

C 


57 

C 


CALL MOVE I 
CALL MOVEE 

Write the Electric field, Potential and Charge Density at 
center point of the system at instance time. 

WRITE ( 92, 41) TIME, EX( 65, 33) ,EY(65,33) 

WRITE (93, 42) TIME, PHI (65, 33) 

WRITE ( 94 , 42 )TIME,RHO( 65, 33 ) 
format (Ix,f8.3,2el2.3) 
format (Ix,f8.3,el2.3) 


I F ( LAC .EQ. 2500 ) THEN 

ITIM = ITST + ( LAC*LWRTE)/5 
CALL WRITEFIL( ITIM) 

LWRTE » LWRTE+1 
LAC = 0 
END IF 

CONTINUE 

CLOSE( 92 ) 

CLOSE (93) 

CLOSE (94) 

STOP 

END 

************************************************************ 
SUBROUTINE INITIAL 

************************************************************ 
REAL ION 

COMMON/PAR/E LE ( 163840, 5), I ON (163840, 5) 

COMMON/ARRAY 1 /PH I (129,65) 

COMMON/EFI ELD/EX (129,65) ,EY( 129 ,65 ) 

COMMON/ARRAY/RHO ( 129,65) , RHOE( 129 , 65 ) , RHOI ( 129 , 6 5 ) 

COMMON/D EL/DVX I , VXIW( 301 ) , FIW( 301 ) , GIW( 301 ) , VDRI , DVXE , VDRE 

> ,VXIO( 301 ) , FIO( 301 ) ,GIO( 301 ) 

> , VXEO( 241 ) , FEO( 241 ) , GEO( 241 ) 

COMMON/VEX/XX ( 1 2 9 ) , YY ( 6 5 ) 

COMMON TIME , NI J 

COMMON NX , NY ,NXl ,NYl , PI , MAIN 

COMMON RATO , RATW , DELT , VELVAR , BXM , BYM , BZM , EYO 
COMMON I SEED , DY , DX ,NIED,NIED1, NSUMO 

OPEN ( UNI T=8 5 , FILE= ' e-v_t0 . dat ' , STATUS = ' NEW ' ) 

OPEN ( UNIT=8 6 , FILE= ' e-p_t0 . dat ' , STATUS= 'NEW' ) 

OPEN ( UNI T=9 5 , FILE= ' i-v_t0 . dat ' , STATUS='NEW' ) 

OPEN ( UNI T=9 6 , F I LE= ' i-p_t0 .dat' , STATUS='NEW' ) 

DO 55 L=1 , NI ED 

WRITE ( 9 5, 19 ) ION (L, 3) , ION ( L , 4 ) , ION ( L , 5 ) 

FORMAT ( IX, 3E12 . 3 ) 

* CONTINUE 

DO 56 L=1 , NI ED 

WRITE (96,17) ION ( L , 1 ) ,ION(L,2) 

FORMAT ( IX, 2F9 . 3 ) 

CONTINUE 

DO 57 L=1 , NI ED 

WRITE (86,17)ELE(L,1) ,ELE(L,2) 

CONTINUE 

DO 58 L= 1 , NI ED 

WRI TE(85,19)ELE(L,3) , ELE ( L , 4 ) , ELE ( L , 5 ) 


the 


****** 

****** 



n O 


CONTINUE 

CLOSE (95) 
CLOSE( 96 ) 
CLOSE ( 85 ) 
CLOSE( 86 ) 
RETURN 
END 


^ 

> ^ HA RAC TER * 1 5 *NAME2 *NAME 3 ^ NAME 4 ^ NAME 5 * NAME 6 *NAME **************** 
REAL ION NAME8 ,NAME9 ,NAME10 ,NAME11 

?omS^o^ E(163840 ' 5) ' ION ( 1 63840,5) 

COMMON/ARRAYl/PHI ( 129 , 65 ) 

CoS/AR^Y/BHo! HI II ) [lloluiVh) RHOK129 65) 

co HMON/DEI ./ovxx , VXXW, 3 0, , , 30 , ; 5 , V DRE 

' / VXEO( 241 ) , FEO( 241 ) GEOf?d1i 

COMMON/VEX/XX (129), Y Y (65) ' 

COMMON TIME , NI J 

COMMON NX, NY, NX1 , NY1 , PI , MAIN 

COMMOM ™2 A**™ ' DELT ' VELVAR ' B XM r BYM , BZM , EYO 
: COMMON I SEED , DY , DX ,NIED,NIED1 , NSUMO 

WRITE ( NAME2 , ' ( ' ' phi t",i4)') ITIM 
OPEN ( UNI T-2 , FI LE=NAME2 , STATUS= ' NEW' ) 

WRITE ( NAME 3 , ' ( • ' rhoi t",i4)') itim 
OPEN( UNIT=3 , FILE=NAME3 , STATUS= 'NEW' ) 

■ WRITE ( NAME4, '( " rhoe t",i4)') ITIM 
. OPEN( UNIT=4 , F I LE=NAME4 , STATUS = ' NEW' ) 

WRITE ( NAME 5 , ' ( ' ' rho t",i4)') ITIM 
OPEN( UNIT=5 , FILE=NAME5 , STATUS= ' NEW' ) 


C 


c 


c 


c 


c 


c 


c 


WRITE ( NAME6, '( "ex t",i4)') ITIM 
OPEN ( UNI T=6 , FI LE=NAME6 , STATUS= ' NEW' ) 

WRITE ( NAME 7 , ' ( " ey t",i4)') ITIM 
OPEN( UNIT-7 , FI LE=NAME7 , STATUS = ' NEW' ) 

WRITE( NAME 8 , ' ( " e-p t",i4)') ITIM 
OPEN( UNIT=8 , F I LE=NAME8 , STATUS= 'NEW' ) 

WRITE ( NAME9, '( "e-v t",14)') ITIM 
OPEN( UNIT=9 , F I LE=NAME9 , STATUS= ' NEW' ) 


WRITE ( NAME1 0 , ' ( " i-p t" 
OPEN ( UNI T=1 0 , FI LE=NAME1 0 


, i4 ) ' ) , ITIM 
, STATUS= ' NEW' ) 


WRITE ( NAME1 1 
OPEN(UNIT=ll 


' ( ' ' i-v_t ' ' , i 4 ) ' ) ITIM 
FILE=NAMEll , STATUS = ' NEW' ) 


WRITE ( NAME2 , ' ( ' 'phi t" ,i4) ' ) 
OPEN ( UNI T=2 , F I LE=NAME2 , STATUS = 


ITIM 
' NEW' ) 


DO 10 J=1 , NY 

WRITE (2,112) ( PHI ( I , J ) 
WRITE ( 3 , 1 1 2 ) ( RHO I ( I , J ) 
WRITE ( 4 , 112 ) ( RHOE ( I , J ) 


1=1 , NX) 

1 = 1 , NX) 
1=1 , NX) 



112 

10 


C 


20 

C 


118 

25 

C 


30 

C 


35 


//'RITE (5,112) ( RHO ( I , J 
WRITE (6,112) ( EX ( I , J ) 
WRITE (7,112)(EY(I,J) 
FORMAT ( IX, 6E12 . 3 ) 
CONTINUE 
CLOSE ( 2 ) 

CLOSE ( 3 ) 

CLOSE ( 4 ) 

CLOSE ( 5 ) 

CLOSE ( 6 ) 

CLOSE ( 7 ) 


I = 1 , NX ) 
1=1 , NX) 

1 = 1 , NX) 


DO 20 L=1 , NIED 
WRITE ( 8 , 1 1 3 ) ELE 
CONTINUE 


L,l) , ELE ( L , 2 ) 


DO 25 L=1 , NIED 

WRITE (10,118) ION ( L , 1 
format ( lx, 2f9 . 3 ) 
CONTINUE 


) , ION ( L , 2 ) 


DO 30 L=1 , NIED 

WRITE (9,114)ELE(L,3) 
CONTINUE 


ELE ( L , 4 ) , ELE ( L , 5 ) 


DO 35 L=1 , NIED 

WRITE ( 11,114) ION ( L , 3 ) 
CONTINUE 


, ION ( L , 4 ) , ION(L, 5 ) 


113 

114 
C 


c 


FORMAT ( IX , 2F9 . 3 ) 
FORMAT ( IX, 3E12 . 3 ) 

CLOSE ( 8 ) 

CLOSE ( 9 ) 

CLOSE ( 10 ) 

CLOSE (11) 

RETURN 

END 


C 

c 


c 

c 

c 

c 


- 

**re1l*i*on***********************************************‘' 
COMMON/PAR/E LE ( 163840, 5), I ON (163840 5) 

COMMON/APRAYl/PHI ( 129,65) UNlI ^ 04U ' b) 

COMMON/EFI ELD/EX ( 129,65) ,EY( 129,65) 

COMMON/ARRAY/RHO (129,65) , RHOE (129 65) RHOK129 API 
COMMON/DEL/DVXI , VXIW ( 3 01 EIw! 301 i “K ! ^ JvXE. VDP 
, VXIO( 301 ) , FIO( 301 ) ,GIO( 301 ) 

, VXEO (241), FEO (241) , GEO (241) 

COMMON/VEX/XX (129) , YY (65) ' 

COMMON TIME , NI J 

COMMON NX , NY , NX1 , NYl , PI , MAIN 

rnMMOM rtc2A RATW ' DELT ' VELVAR , BXM , B YM , BZM , EYO 
COMMON I SEED , DY , DX , NIED , NI EDI , NSUMO 

COMMON COKX ( 129 ) , COKY( 65) , SM( 129 , 65 ) 


PI=3 .1415926535898 
Normalized Magnetic field 
BM=0 . 3 

Angle of the magnetic field to 
THETA=0 .05 

Magnetic field in Z-component 
BZM=BM*COS ( THETA ) 


the x-axis 


* * * 
★ * * 



n n n 


c 

C 

c 


c 


1 

2 

c 

10 

c 

20 


30 


Magentic field in 
BYM=0 . 


Y-component 


Magnetic field in 
BXM=BM*SIN( THETA) 
Normalized Water I 
VELVAR =0.2 
EYO = 0. 

DX = 2. 

DY = 2 . 

RATW = 1.0/(2304.) 
RATO = 1 . 0/( 2048 . ) 


X-component 
ons Drift veloci 


ty 


DO 1 J=1 , NY 

YY( J) =DY* FLOAT ( J-l ) 
CONTINUE 
DO 2 1=1, NX 

XX( I ) = DX* FLOAT ( 1-1 ) 
CONTINUE 


DO 10 1=1, NX1 
COKX(I) = ( 
CONTINUE 


2 . *PI*I/XX( NX) 


) **2 


DO 20 1=1, NY1 
COKY(I) = ( 
CONTINUE 


2 . *PI * I/YY ( NY ) 


) **2 


Smoothing Coefficient 


DO 30 J=1 , NYl 
DO 30 1=1, NX1 
SM(I,J) = (NX1*SIN(PI* 
CONTINUE 


I/NX1 ) *NY1 


* S I N ( P I * J/NY 1 )/(PI*I*PI*J)) 


RETURN 

END 


C 

C 


C 


10 

c 


“”»^»r;,*;;,*^****************“***»**»»*»*»**» 

**reIl*i*on**********************‘**************************** 

COMMON/EFIELD/EX ( 129, 65), EY (129 65) 

COMMON/ARRA Y/RHO (12 9,65) , RHOE (129 65) RHOK129 fi?) 

, VXEO( 241 ) , FEO( 241 ) GEO( 241 ) 

COMMON/VEX/XX (129) , Y Y (65) 

COMMON TIME , NI J 

COMMON NX , NY , NX1 , NYl , PI , MAIN 

COMMON RATO, RATW, DELT, VELVAR, BXM,BYM BZM EYO 
COMMON I SEED , DY , DX , NI ED , NI EDI , NSUMO ' ' 


DVXI=0 .00104 
VDRI = VELVAR 
VXIW(l) = 0. 

FIW(l) = 0. 

GIW(l) = 0. 

DO 10 I = 2,301 

VXIW(I) = FLOAT ( 
GIW(I) =EXP( - ( 57 
FIW(I) = DVXI * 0 . 
CONTINUE 


1-1 ) *DVXI + VXIW ( 1 ) 

6 . * 4 . * ( VXIW ( I ) -VDRI ) * * 2 ) /2 ) 
5* ( GIW ( 1-1 )+GIW( I ) ) +FIW ( 1-1 ) 


DO 20 J = 
F IW ( J 


1 ,301 

= FIW( J )/FIW( 301 ) 


**2 

★ 

★ 



o n 


20 


10 


CONTINUE 

RETURN 

END 

aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa*aaaaaaaaaaaaaaaaaaaaaaaaa 

SUBROUTINE DISTELE 

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
REAL ION 

COMMON/PAR/ELE( 163840 , 5 ) , ION( 16 38 40 , 5 ) 

COMMON/ARRAYl/PHI (129,65) 

COMMON/E FI ELD/EX ( 129 , 65 ) , EY( 129 , 65 ) 

rnMUnM^f^ Y/RH0( 129,65 } ' RH0E ( 129 ' 65 ) , RHOI ( 129 , 65 ) 
COMMON/DEL/DVXI , VXIW( 301) ,FIW( 301 ) ,GIW(301) , VDRI , DVXE VDRE 

> , VX IO( 301) , F I O ( 3 0 1 ) ,GIO( 301) 

> , VXEO (241),FEO(241), GEO (241) 

COMMON/VEX/XX ( 1 2 9 ) , Y Y ( 6 5 ) 

COMMON TIME , NI J 

COMMON NX , NY , NXl , NYl , PI , MAIN 

COMMON RATO , RATW , DELT , VELVAR , BXM , BYM , BZM , EYO 
COMMON I SEED , DY , DX ,NIED,NIED1, NSUMO 

DVXE- 0.05 
VDRE - 0. 

VXEO ( 1 ) = -6. 

FEO ( 1 ) - 0. 

GEO ( 1 ) . 0. 

DO 10 I = 2,241 

VXEO ( I ) = FLOAT( 1-1 ) *DVXE +~VXEO(l) 

GEO ( I ) «EXP(-( ( VXEO ( I ) -VDRE ) * * 2 ) /2 ) 

CONTINUE^ = DVXE *°- 5 * (GE °( i - 1 ) +ge:o ( i ) )+FEO(I-l) 


20 


DO 20 J - 1,241 

FEO(J) = FEO ( J ) /FEO (241) 
CONTINUE 
RETURN 
END 


C 

C 


10 


AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
SUBROUTINE DISTIONO 

H**********************************************************. 

REAL ION 

COMMON/PAR/ELE ( 163840,5) , I ON (163840, 5) 

COMMON/ARRAYl/PHI (129,65) 

COMMON/EFI ELD/EX ( 129, 65), EY (129, 65) 

COMMON/ARRAY/RHO (129,65), RHOE (129,65), RHOI (129,65) 
COMMON/DEL/DVXI VXIW( 301 ) , FIW( 301 ) ,GIW< 301 ) , VDRI , DVXE , VDRE 

> , VX IO( 301) ,FIO( 301 ) ,GIO( 301) 

> , VXEO( 241 ) , FEO( 241 ) ,GEO( 241 ) 

COMMON/VEX/XX ( 1 2 9 ), YY ( 6 5 ) 

COMMON TIME , NI J 

COMMON NX, NY, NXl, NYl, PI, MAIN 

COMMON RATO , RATW , DELT , VELVAR , BXM , BYM , BZM , EYO 
COMMON ISEED,DY,DX,NIED,NIED1, NSUMO 

DVXI=0 .0011 
VDRI = 0. 

VXIO(l) = -0.1325 
FIO(l) = 0. 

GIO(l) = 0. 

DO 10 I = 2,301 

VXIO(I) = FLOAT ( I -1 ) *DVXI + VXIO(l) 

VXIO(I) = FLOAT ( 1-1 ) *DVXI 
GIO ( I ) =EXP ( - ( 5 1 2 . * 4 . * ( VX IO(I) -VDRI ) **2 )/2 ) 

CONTINUE^ = DVXI * 0 - 5 * (GIO ( i - 1 ) +gi °( i ) )+FIO(I- l) 


* * * 
***c 


* * 
* ★ 



non 


DO 20 J = 1,301 

FIO(J) = FIO( J)/FIO( 301 ) 
CONTINUE 
RETURN 
END 


*‘“M^w“Kir*****“*************“**“*“*“*“****»** 

**R£^L*jqjj** ******************************** * ******************* 

COMMON/PAR/ELE( 163840, 5) ,I0N(16 3840 5) 

COMMON/ARRAYl/PHI( 129,65 ) U ^ 1Dja<1U ' b > 

COMMON/EFIELD/EX ( 129 , 65 ) ,EY(129 65) 

COMMON/ARRAY/RHO (12 9,65) RHOE ( 1 2 Q fi c: \ BHnT , no c c , 
COMMON/DEL/DVXI j j , FIW( 301 j ^GIWUOl ^VORI f DVXE , VDRE 

__ > VXEO ( 2 41 ) , FEO ( 241 ) , GEO( 2 4 1 ) 

COMMON/VEX/XX (129) , Y Y (65) 

COMMON TIME , NI J 

COMMON NX, NY, NXl ,NY1 , PI , MAIN 

COMMOM ' DELT ' VELVAR ' BXM , BYM , BZM , EYO 

COMMON I SEED , DY , DX, NIED , NIED1 , NSUMO 

PI= 3.1415926535898 

TPI = 2.*PI 

Distribute 20 electrons per cell in the system 

NSUM = -3 
DO 13 J=1 , NY1 
DO 14 1=1, NXl 
DO 12 M-1,3 

NSUM = NSUM + 3 
DO 11 L=1 , 3 

ELE(L+NSUM,1)= 0 . 5*L + DX* ( I— 1 ) 

ELE ( L+NSUM , 2 ) = 0.5*M + DY*(J-1) 

RND = RAN ( I SEED ) ’ 

DO 40 K=1 ,241 
DIFFl=FEO( K) -RND 
I F ( DI FF1 .GE. 0.0) GOTO 20 
CONTINUE 

ELE ( L+NSUM , 3 ) = VXEO(K) 

= 1 . 4 J. 4 * SQRT ( - LOG ( 1-RAN ( I SEED ) ) ) 

TPAR = TPI*RAN( ISEED) 

ELE( L+NSUM, 4 ) = VPAR* S IN ( TPAR ) 

ELE ( L+NSUM , 5 ) = VPAR*COS ( TPAR ) 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

DO 133 J=1 , NY1 
DO 144 1=1, NXl 
DO 122 M-1,3 

NSUM = NSUM + 3 
DO 111 L-1,3 

ELE( L+NSUM, 1 ) = 0.5*L + DX*(I-1) + 0 25 
ELE ( L+NSUM , 2 ) = 0.5*M + DY* ( J-l ) - 0 25 
RND = RAN ( ISEED ) 

DO 440 K=1 ,241 
DIFFl=FEO(K)-RND 
I F ( DI FF1 .GE. 0.0) GOTO 220 
CONTINUE 

ELE ( L+NSUM , 3 ) = VXEO(K) 

VPAR = 1.414* SQRT ( -LOG ( 1-RAN ( ISEED ) ) ) 



Ill 

122 

144 

133 

c 


484 

282 


81 

82 

84 

83 

C 


TPAR = TPI*RAN( ISEED) 

ELE ( L+NSUM , 4 ) = VPAR*SIN( TPAR ) 
ELE ( L+NSUM , 5 ) = VPAR*COS ( TPAR ) 
CONTINUE 
CONTINUE 
CONTINUE 
CONTINUE 


1 ) ) 

. EQ. 


1 ) .AND. (M . EQ. 
■EQ. 3 ) . AND . ( M 
+ 1 

0 . 5*L + DX* ( 1-1 ) 
,5*M + DY* ( J-l ) 


0 


NSUM = NSUM+3 
DO 83 J=1 , NY1 
DO 84 1=1, NX1 
DO 82 M-1,3 
DO 81 L=1 , 3 
IF( ( ( L . EQ. 

.OR. ((L 
NSUM = NSUM 
ELE (NSUM, 1 )= 

ELE ( NSUM , 2 ) = 

RND = RAN (I SEED) 

DO 484 K=1 ,241 
DIFFl=FEO(K)-RND 
I F ( DI FF1 .GE. 0.0) GOTO 282 
CONTINUE 

ELE (NSUM, 3) = VXEO(K) 

VPAR = 1.41 4 *SQRT ( -LOG ( 1-RAN ( I SEED 
TPAR = TPI*RAN( ISEED) 

ELE(NSUM,4) = VPAR*SIN( TPAR ) 

ELE (NSUM, 5) = VPAR* COS ( TPAR ) 

END IF 
CONTINUE 
CONTINUE 
CONTINUE 
CONTINUE 


3 ) ) ) THEN 

- 0.25 
+ 0.25 


) ) ) 


C 

c 

c 


409 

209 


511 

312 
314 

313 
C 


PRINT*, 'Total Number of Electrons in the system 
Distribute 10 Oxygen Ions per cell in the sysytem 


NSUMO = 0 
DO 313 J=1 , NY1 
DO 314 1=1, NX1 
DO 312 M=l,3 
DO 511 L=1 , 3 

NSUMO = NSUMO + 1 
ION ( NSUMO , 1 ) = 0 . 5 * L 
I ON ( NSUMO , 2 ) = 0 . 5 *M 
RND = RAN ( ISEED ) 

DO 409 K=1 ,301 
DIFFl=FIO( K) -RND 


+ DX* ( I — 1 ) 
+ DY* ( J— 1 ) 


+ 


I F ( DI FFl .GE. 0.0) GOTO 209 
CONTINUE 

ION ( NSUMO , 3 ) = VXIO(K) 

VPAR = 0 . 03125*SQRT( -LOG ( 1-RAN ( ISEED 
TPAR = TPI*RAN( ISEED) 

ION ( NSUMO , 4 ) = VPAR*SIN( TPAR) 

ION ( NSUMO , 5 ) = VPAR*COS ( TPAR ) 
CONTINUE 
CONTINUE 
CONTINUE 
CONTINUE 


0.25 

0.25 


) ) 


DO 73 J=1 , NY1 
DO 74 1=1, NX1 
DO 72 M-1,3 
DO 71 L-1,3 
IF((L .EQ. 


1) .AND. (M .EQ. 1 ) ) THEN 



c 


479 

279 

C 

C 

C 

C 


71 

72 
74 

73 
C 

c 

c 

c 


461 

261 


711 

712 
714 

713 
C 


661 

264 


NSUMO = NSUMO + 1 

ION ( NSUMO , 1 ) = 0 . 5 * L + DX*(I-1) - 0.25 
I ON ( NSUMO , 2 ) = 0 . 5 *M + DY * ( J-l ) + 0.25 

RND = RAN (I SEED) 

DO 479 K-1,301 
DIFFl=FIO( K ) -RND 
I F ( DI FFl . GE . 0.0) GOTO 279 
CONTINUE 

ION (NSUMO, 3) = VXIO(K) 

SQRT(IOO) =10. 

1/10. = 0.1 
1.414 * 0.1 = 0.1414 

VPAR = 0 . 0666667 *SQRT( -LOG( 1-RAN( I SEED ) ) ) 

VPAR = 0.0312 5 * SQRT ( -LOG ( 1-RAN ( I SEED ) ) ) 

TPAR = TPI*RAN( ISEED) 

ION ( NSUMO , 4 ) = VPAR* S IN ( TPAR) 

ION ( NSUMO , 5 ) = VPAR*COS ( TPAR ) 

ENDIF 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

PRINT* , ' Total Number of Oxygen Ions in the system= ' , NSUMO 

Distribute 10 Water Ions per cell in the system 

NSUM = NSUMO 
DO 713 J=1 , NYl 
DO 714 1=1, NX1 
DO 712 M=1 , 3 
DO 711 L=1 , 3 

NSUM = NSUM + 1 

ION ( NSUM , 1 ) = 0 . 5 *L + DX* ( 1-1 ) - 0.125 
ION ( NSUM , 2 ) = 0 . 5 *M + DY*(J-1) + 0.125 
RND = RAN( ISEED) 

DO 461 K-1,301 
DIFFl=FIW(K) -RND 
I F ( DI FFl .GE. 0.0) GOTO 261 
CONTINUE 

ION ( NSUM , 3 ) = VXIW(K) 

VPAR = 0.029 5 * SQRT ( -LOG ( 1-RAN ( ISEED ) ) ) 

TPAR = TPI *RAN( ISEED ) 

ION ( NSUM , 4 ) = VPAR*SIN( TPAR) 

ION (NSUM, 5) = VPAR*COS ( TPAR ) 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

DO 613 J-l, NYl 
DO 614 1=1, NX1 
DO 612 M= 1 , 3 
DO 611 L=1 , 3 

I F ( ( L . EQ. 1 ) . AND . ( M . EQ . 1 ) ) THEN 

NSUM = NSUM + 1 

ION ( NSUM , 1 ) = 0 . 5*L + DX*(I-1) + 0.125 
I ON ( NSUM , 2 ) = 0 . 5 *M + DY*(J-1) - 0.125 
RND = RAN ( ISEED) 

DO 661 K-1,301 
DIFFl=FIW( K) -RND 
I F ( DI FFl .GE. 0.0) GOTO 264 
CONTINUE 

ION ( NSUM , 3 ) = VXIW(K) 

VPAR = 0.0295 * SQRT ( -LOG ( 1-RAN ( I SEED ) ) ) 



611 

612 

614 

613 

C 


C ** 

c ** 


c 


5 

6 
C 


C 

C 


C 


1 

c 


TPAR = TP I *RAN ( ISEED ) 

ION ( NSUM , 4 ) = VPAR*SIN( TPAR) 

ION ( NSUM , 5 ) = VPAR*COS( TPAR) 

ENDIF 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

PRINT* , ' Total Number of Ions in the system® NSUM 

NIED = NSUM 

NIED1 = NIED -1 

RETURN 

END 


*************************************************************** 
SUBROUTINE DENS 

*************************************************************^^ 
INTEGER S IGR , S IGT 
REAL ION 

COMMON/PAR/ELE( 16 38 4 0, 5) , ION( 16 3840 , 5 ) 

COMMON/ARRAY1/PHI (129,65) 

COMMON/EFI ELD/EX ( 129,65) ,EY(129,65) 

COMMON/ARRAY/RHO (129,65), RHOE ( 1 2 9 , 6 5 ) , RHOI ( 1 2 9 , 6 5 ) 

COMMON/D EL/DVX I , VXIW ( 301 ) ,FIW( 301 ) ,GIW( 301) , VDRI , DVXE , VDRE 

> ,VXIO( 301 ) ,FIO( 301 ) ,GIO( 301 ) 

> , VXEO( 241 ) , FEO( 241 ) ,GEO( 241 ) 

COMMON/VEX/XX ( 1 2 9 ) , YY ( 6 5 ) 

COMMON TIME , NI J 

COMMON NX , NY , NX1 , NYl , PI , MAIN 

COMMON RATO , RATW , DELT , VELVAR , BXM , B YM , BZM , EYO 
COMMON ISEED, DY,DX, NIED, NIEDl ,NSUMO 

DO 6 J=1 , NY 

DO 5 1=1, NX 
RHOE ( I , J ) =0 . 

RHOI( I, J)=0. 

RHO ( I , J ) =0 . 

CONTINUE 

CONTINUE 

DO 1 L= 1 , NIED 

NI I =INT ( ELE ( L , 2 ) /DY ) +1 
I F ( NI I .EQ. NY+1 ) NI 1=1 
NJ=INT(ELE(L,1)/DX)+1 

S IGR=+1 
S IGT=+1 

DELA=ABS ( ELE ( L , 1 ) -XX ( N J ) ) 

DELB=ABS ( ELE ( L , 2 ) - YY (Nil) ) 

Al=DELA*DELB 
A2=DX*DELB-A1 
A3=DY*DELA-A1 
TOTA=DX*DY 
A4=TOTA- ( A1+A2+A3 ) 


RHOE( NJ+SIGR, NII+SIGT ) =RHOE ( NJ+SIGR , NII+SIGT ) +A1/TOTA 

RHOE ( NJ , NI I+S IGT ) =RHOE( NJ , NII+SIGT ) +A2/TOTA 

RHOE ( NJ+SIGR , NI I ) =RHOE ( N J+S IGR , NI I ) +A3/TOTA 

RHOE ( NJ , NI I ) =RHOE ( N J , NI I ) +A4/TOTA 

CONTINUE 

DO 2 L=1 , NIED 

NI I = INT ( I ON ( L , 2 ) /DY ) +1 
I F ( Nil .EQ. NY+1 ) NI 1=1 


C 



n n 


2 

C 


100 

c 


108 

C 


11 

10 

c 


8 

4 

C 


C 

c 


NJ=INT ( ION ( L , 1 ) /DX ) +1 

SIGR=+1 
S IGT=+1 

DELA=ABS ( ION ( L , 1 ) -XX ( NJ ) ) 

DELB=ABS ( ION ( L , 2 ) -YY (Nil) ) 

Al =DELA*DELB 

A2=DX*DELB-Al 

A3=DY*DELA-Al 

TOTA=DX*DY 

A4=TOTA-( A1+A2+A3 ) 


RHOI ( NJ+SIGR, NII+SIGT ) =RHOI ( NJ+SIGR , NI I+SIGT ) +A1/TOTA 

RHOI ( NJ , NI I+SIGT ) =RHOI ( NJ , NII+SIGT ) +A2/TOTA 

RHOI (NJ+SIGR, Nil )=RHOI ( NJ+SIGR, Nil ) +A3/TOTA 

RHOI ( NJ , NI I ) =RHOI ( NJ , NI I ) +A4/TOTA 

CONTINUE 

DO 100 1=1, NX 

RHOI (1,1) =RHOI (1,1) +RHOI ( I , NY ) 

RHOI (I, NY) -RHOI (1,1) 

RHOE ( I , 1 )=RHOE( I , 1 ) +RHOE ( I , NY ) 

RHOE ( I , NY ) =RHOE (1,1) 

CONTINUE 

DO 108 J=1 , NY 

RHOI ( 1 , J ) =RHOI ( 1 , J ) +RHOI ( NX , J ) 

RHOI ( NX , J ) =RHOI ( 1 , J ) 

RHOE ( 1 , J ) =RHOE ( 1 , J ) +RHOE ( NX , J ) 

RHOE ( NX , J ) =RHOE ( 1 , J ) 

CONTINUE 

DO 10 J=1 , NY 

DO 11 1=1, NX 

RHOI ( I , J ) = RHOI ( I , J)/NIJ 
RHOE ( I , J ) = RHOE ( I , J ) /NI J . 

CONTINUE 

CONTINUE 

DO 4 J=1 , NY 
DO 8 1=1, NX 

RHO ( I , J ) = RHOI ( I , J )-RHOE( I , J) 

CONTINUE 

CONTINUE 

RETURN 

END 

*************************************************************** 
SUBROUTINE POTENT 

*************************************************************** 
INTEGER NRA , NCA , LDA , LDCOEF 
COMPLEX RHOCOMP( 129,65) , PHI COMP ( 1 29 , 6 5 ) 

REAL ION 

COMMON/PAR/ELE( 16 3840 , 5), ION (163840, 5) 

COMMON/ARRAY 1 /PH I ( 129 , 65 ) 

COMMON/EFI ELD/EX ( 129,65) ,EY( 129,65) 

COMMON/ARRAY/RHO (129,65), RHOE (129,65), RHO 1(129,65) 

COMMON/D EL/DVX I , VXIW( 301 ) , FIW( 301 ) , GIW( 301 ) , VDRI , DVXE , VDRE 

> ,VXIO( 301) ,FIO( 301) ,GIO( 301) 

> , VXEO( 241 ) , FEO( 241 ) , GEO ( 241 ) 

COMMON/VEX/XX ( 1 2 9 ) , YY ( 6 5 ) 

COMMON TIME , NI J 

COMMON NX , NY , NXl , NYl , PI , MAIN 

COMMON RATO , RATW , DELT , VELVAR , BXM , BYM , BZM , EYO 
COMMON I SEED , DY , DX , NI ED , NI EDI ,NSUMO 



UUUO fMU ^rrou -H U U U 


COMMON COKX (129) ,COKY(65) ,SM( 129,65) 

DIMENSION NN( 2 ) ,DATA( 16384 ) 

C DATA ( 2 *NXl *NYl ) 

PI = 3.1415926535898 
TPI = 2 . *PI 
NN ( 1 ) = NXl 
NN ( 2 ) = NYl 

C SKIP 1 = 1 ( LENGHT=0 ) BECAUSE IT SHOULD START FROM LENGTH= 1 

L= 2 *NXl 
DO 10 J=2 , NY 
DO 10 1=2, NX 

DATA( 2* (1-2) + 1 + L * ( J-2 ) ) = RHO(I,J) 

DATA( 2 * ( 1-1 ) + L * (J-2) ) = 0. 

10 CONTINUE 

Using two dimensioal FFT .( forward ) 

CALL FOURN ( DATA, NN , NDIM , IS IGN ) 

CALL FOURN ( DATA, NN, 2 , 1 ) 

DO 20 J=2 , NY 

DO 20 1 = 2, NX '■ 

REALRHOT = DATA ( 2* (1-2) +1 + L * (J-2) ) 

AIMAGRHOT= DATA ( 2*(I-1) + L * (J-2) ) 

RHOCOMP ( I , J ) = REALRHOT* ( 1 ., 0 . ) + AIMAGRHOT* ( 0 . , 1 . ) 

0 CONTINUE 

DO 30 J=2 , NY 
DO 40 1=2, NX 

PHI COMP ( I , J ) =RHOCOMP ( I , J ) * SM ( 1-1 , J-l )/( COKX( 1-1 ) +COKY( J-l ) ) 
0 CONTINUE 

0 CONTINUE 

DO 100 J=2 , NY 
DO 100 1=2, NX 

DATA ( 2* (1-2) + 1 + L * (J-2) ) = REAL ( PHICOMP ( I , J ) ) 

DATA ( 2* ( 1-1 ) + L * (J-2) ) = AIMAG ( PHICOMP ( I , J ) ) 

00 CONTINUE 

Using two dimensional FFT .( backward ) 

CALL FOURN( DATA, NN, NDIM, ISIGN) 

CALL FOURN ( DATA , NN ,2,-1) 

DO 200 J=2 , NY 
DO 200 1=2, NX 

REALPHIT = DATA ( 2* (1-2) +1 + L * (J-2) ) 

AIMAGPHIT= DATA( 2*(I-1) + L * (J-2) ) 

PHICOMP ( I , J ) = REALPHIT* ( 1 ., 0 . ) + AIMAGPHIT* ( 0 . , 1 . ) 

200 CONTINUE 

C 

DO 50 1=2, NX 
DO 60 J=2 , NY 

C PH I ( I , J ) = REAL ( PHICOMP ( I , J ) ) /( XX ( NX ) * YY ( NY ) ) 

PHI ( I , J ) = REAL ( PHICOMP ( I , J ) ) /( NXl * NYl ) 

60 CONTINUE 

50 CONTINUE 

C 

DO 80 J= 1 , NY 

PHI ( 1 , J ) = PHI ( NX, J ) 

80 CONTINUE 

DO 90 1=1, NX 

PHI (1,1) = PHI ( I , NY ) 

90 CONTINUE 

C 

RETURN 

END 



c 


c 


*********************** 
SUBROUTINE FIELD 


********** ******************** 


* * 

* * 


******** 

******* 


C 


5 

4 

C 


6 

7 

C 


8 


C 


c 


9 


C 


REAL ION 

COMMON/PAR/E LE (163840, 5), I ON (163840,5) 

COMMON/ARRAY1/PHI (129,65) 

COMMON/EFI ELD/EX ( 129,65) , EY( 129 , 65 ) 

rnMMOM^np w Y/RH ° ( 129 ' 6 5 ) > RH°E ( i 29 , 6 5 ) , RHOI (129,65) 

COMMON/D EL/DVX I , VXIW( 301 ) , FIW( 3 01 ) ,GIW< 3 01 ) , VDRI , DVXE , VDRE 
, VXIO( 301), FIO( 301 ) ,GI 0(301) 

, VXEO (241) , FEO (241) ,GEO(241) 

COMMON/VEX/XX ( 1 2 9 ) , YY ( 6 5 ) 

COMMON TIME , NI J 

COMMON NX , NY , NXl , NYl , PI , MAIN 

T ° ' ' DELT ' VELVAR , BXM , B YM , BZM , EYO 

COMMON I SEED , DY , DX , NIED , NIEDl , NSUMO 

DO 4 J=2 , NYl 
DO 5 1=1, NX 

CONTINUE J } = ~ ( PHI ( 1 ' J+1 } _PHI ( 1 ' J_1 > ,/( 2 • °*DY) 

CONTINUE 


DO 7 J=1 , NY 
DO 6 1=2, NXl 

CONTINUE 1 ' J)=_(PHI( 1 + 1 ' J)_PHI( I_1 ' J) ,/( 2 -°* DX) 
CONTINUE 


DO 8 J=1 , NY 

EX (1,J)=— (PHI(2,J)-PHI(1 
EX ( NX , J ) =— ( PHI ( NX , J ) -PHI 
CONTINUE 


r J ) ) /DX 
(NXl, J) ) /DX 


DO 9 1=1, NX 

EY( I , 1 )=-( PHI ( I , 2 ) —PHI ( I 
EY ( I , NY ) =EY (1,1) 

EY( I ,NY)=-( PHI ( I , NY) -PHI 
CONTINUE 


, 1 ) ) /DY 
( I , NYl ) ) /DY 


RETURN 

END 


C 

c 


c 


******************************^^^^^^^^^^^ 
SUBROUTINE MOVEI 

* t III ********************************************************* 

INTEGER S IGT , S IGR 
REAL ION 

COMMON/PAR/E LE (163840,5) , I ON (163840, 5) 

COMMON/ARRAY 1 /PH I (129,65) 

COMMON/EFI ELD/EX ( 129,65) , EY( 129 , 65 ) 

COMMON/ARRAY/RHO (129,65) , RHOE (129,65), RHOI (129 65) 

COMMON/D EL/DVX I , VXIW ( 301) , FIW( 3 d ) ,GIW(301) , VDRI , DVXE , VDRE 

> , VX 10(301) ,FIO( 301) , G I O ( 3 0 1 ) 

> , VXEO (241) , FEO ( 2 4 1 ) , GEO ( 2 4 1 ) 

COMMON/VEX/XX ( 1 2 9 ), YY ( 6 5 ) 

COMMON TIME , NI J 

COMMON NX, NY, NXl, NYl, PI , MAIN 

COMMON RATO , RATW , DELT , VELVAR , BXM , BYM , BZM , EYO 
COMMON I SEED , DY , DX , NIED , NIEDl , NSUMO 


DO 1 L=1 , NIED 

I F ( I ON ( L , 1 ) .GT. XX ( NX ) ) GO TO 1 
I F ( I ON ( L , 1 ) .LT. XX ( 1 ) ) GO TO 1 
NI 1= INT ( ION ( L , 2 ) /DY ) +1 



N J= I NT ( ( I0N(L,1) )/DX)+l 

S IGR=+1 . 

SIGT=+1 . 

DELA=ABS ( ION ( L , 1 ) -XX ( NJ ) ) 

DELB=ABS( ION ( L , 2 )-YY(NII ) ) 

Al=DELA*DELB 
A2=DX*DELB-Al 
A3=DY*DELA-A1 
A4=DX*DY-( A1+A2+A3 ) 

TOTA=DX*DY ;< 

EXPF=» ( Al *EX ( NJ+SIGR , NI I+S IGT ) +A2 *EX ( NJ , NI I +S I GT ) + 

> A3*EX( NJ+SIGR, Nil ) +A4 *EX( NJ , NI I ) )/TOTA 

EYPF=( Al*EY( NJ+SIGR, NI I+S IGT) +A2*EY(NJ,NI I+S IGT) + 

> A3*EY( NJ+SIGR, Nil ) +A4 * EY ( NJ , NI I ) )/TOTA 
C 

I F ( L .LE. NSUMO) THEN 

ION ( L , 3 ) =ION ( L , 3 ) + ( EXP F+ I ON ( L , 4 ) *BZM-ION ( L , 5 ) *BYM ) *DELT*RATO 
ION( L , 4 ) =ION ( L , 4 ) + ( EYO+EYPF- I ON ( L , 3 ) *BZM+ION( L, 5 ) *BXM ) *DELT*RATO 
ELSE 

I ON ( L , 3 ) = I ON ( L , 3 ) + ( EXPF+ION( L, 4 ) *BZM-ION( L , 5 ) *BYM ) *DELT*RATW 
ION ( L , 4 )=ION( L, 4 )+( EYO+EYPF- ION ( L, 3 ) *BZM+ION( L, 5 ) *BXM ) *DELT*RATW 
END IF 

C 

XAX = ION( L , 1 ) 

YAY = ION ( L , 2 ) 

C 

I ON ( L , 1 ) = I ON ( L , 1 ) + I ON ( L , 3 ) *DELT 
ION ( L , 2 ) = I ON ( L , 2 ) + I ON ( L , 4 ) *DELT 
I F ( L .LE. NSUMO) THEN 

I ON ( L , 5 ) = I ON ( L , 5 ) + ( ION ( L , 1 ) -XAX ) *BYM*RATO 

> - ( ION ( L , 2 ) -YAY ) *BXM*RATO 
ELSE 

I ON ( L , 5 ) = I ON ( L , 5 ) + ( ION ( L , 1 ) -XAX ) *BYM*RATW 

> - ( ION ( L , 2 ) -YAY ) *BXM*RATW 
END IF 

1 CONTINUE 

C 

DO 111 L=1 , NI ED 

I F ( ION ( L , 1 ) .GT. XX (NX) ) THEN 
ION ( L , 1 ) = ION ( L , 1 ) - XX (NX) 

ELSE 

IF( I ON ( L , 1 ) .LT. 0.0 ) THEN 

ION ( L , 1 ) = XX ( NX ) + I ON ( L , 1 ) 

ENDIF 
END IF 

111 CONTINUE 
C 

DO 112 L=1 , NI ED 

I F ( ION ( L , 2 ) .GT. YY(NY) ) THEN 
I ON ( L , 2 ) = I ON ( L , 2 ) - YY(NY) 

ELSE 

I F ( ION ( L , 2 ) .LT. 0.0 ) THEN 

ION ( L , 2 ) = YY(NY) + I ON ( L , 2 ) 

ENDIF 

ENDIF 

112 CONTINUE 

C 

RETURN 

END 

C *************************************************************** 

SUBROUTINE MOVEE 

C *************************************************************** 

INTEGER SIGT , SIGR 
REAL ION 

COMMON/PAR/ELE (163840,5) ,ION( 163840,5) 

COMMON/ARRAY 1 /PH I (129,65) 



COMMON/EF I ELD/EX (12 9,65) ,EY( 129 ,65) 

COMMON/ARRAY/RHO ( 129,65) ,RHOE( 129,65) , RHOI ( 129 , 65 ) 
COMMON/DEL/DVXI , VXIW( 301 ) , FIW( 301 ) ,GIW( 301 ) , VDRI , DVXE , VDRE 

> ,VXIO( 301 ) , FIO( 301 ) ,GIO( 301 ) 

> , VXEO( 241) , FEO( 241 ) , GEO( 241 ) 

COMMON/VEX/XX ( 129 ) , YY( 65 ) 

COMMON TIME , NI J 

COMMON NX , NY , NXl , NYl , P I , MAIN 

COMMON RATO , RATW , DELT , VELVAR , BXM , BYM , BZM , EYO 
COMMON I SEED , DY , DX , NI ED , NIEDl , NSUMO 
C 

DO 1 L=1 , NIED 

I F ( ELE ( L , 1 ) .GT. XX(NX)) GO TO 1 
I F ( ELE ( L , 1 ) .LT. XX(1)) GO TO 1 
NI 1= INT ( ELE ( L , 2 ) /DY ) +1 
NJ=INT( ( ELE ( L , 1 ) )/DX)+l 
S IGR=+1 . 

SIGT=+1 . 

DELA=ABS ( ELE ( L , 1 ) -XX ( NJ ) ) 

DELB=ABS ( ELE ( L , 2 ) -YY (Nil) ) 

Al=DELA*DELB 
A2=DX*DELB-A1 
A3=DY*DELA-Al 
A4-DX*DY-( A1+A2+A3 ) 

TOTA=DX*DY 

EXPF = ( Al * EX ( N J + S IGR , NI I + S I GT ) +A2 *EX (NJ ,NII + SI GT ) + 

> A3* EX ( N J+S I GR , NI I ) +A4 *EX(NJ,NII ) J/TOTA 

EYPF* ( Al *EY ( N J+S IGR , NI I+S IGT ) +A2*EY ( NJ , NI I+S IGT ) + 

> A3 *EY ( NJ+S IGR , NI I ) +A4 *EY ( NJ , NI I ) )/TOTA 
C 

ELE (L,3)=ELE(L,3)-( EXP F+ ELE ( L , 4 ) *BZM-ELE(L, 5 ) *BYM ) *DELT 
ELE( L, 4 )=ELE( L, 4 )-( EYO+EYPF-ELE ( L, 3 ) *BZM+ELE ( L , 5 ) *BXM ) *DELT 
C 

XAX = ELE ( L , 1 ) 

YAY = ELE ( L , 2 ) 

ELE( L, 1 )=ELE( L, 1 )+ELE( L, 3 ) *DELT 
ELE ( L , 2 )=ELE( L , 2 ) +ELE ( L, 4 ) *DELT 

ELE( L, 5 )=ELE( L, 5 )-( ( ELE ( L , 1 ) -XAX ) *BYM- ( ELE ( L , 2 ) - YAY ) *BXM ) 

1 CONTINUE 

C 

DO 111 L=1 , NI ED 

I F ( ELE ( L , 1 ) .GT. XX ( NX ) ) THEN 

ELE ( L , 1 ) = ELE ( L , 1 ) - XX(NX) 

ELSE 

I F ( ELE ( L , 1 ) .LT. 0.0 ) THEN 

ELE ( L , 1 ) = XX ( NX ) + ELE ( L , 1 ) 

ENDIF 

ENDIF 

111 CONTINUE 
C 

DO 112 L=1 , NI ED 

I F ( ELE ( L , 2 ) .GT. YY(NY) ) THEN 
ELE ( L , 2 ) = ELE ( L , 2 ) - YY(NY) 

ELSE 

I F ( ELE ( L , 2 ) .LT. 0.0 ) THEN 

ELE ( L , 2 ) = YY(NY) + ELE ( L , 2 ) 

ENDIF 

ENDIF 

112 CONTINUE 
C 

RETURN 

END 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE EXCHANGE 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 



REAL I ONPX , I ONP Y 

REAL ION 

COMMON/PAR/ELE (163840, 5), I ON (16 3840, 5) 

COMMON/ARRAYl/PHI (129,65) 

COMMON/EFI ELD/EX ( 129 ,65),EY(129,65) 

COMMON/ARRAY/RHO (129,65), RHOE (129,65), RHOI (129,65) 

COMMON/D EL/DVX I , VXIW( 301 ) , FIW( 301 ) ,GIW( 301 ) , VDRI , DVXE , VDRE 

> , VXIO( 301 ) , FIO( 301 ) ,GIO( 301 ) 

> , VXEO( 241 ) , FEO( 241 ) ,GEO( 241 ) 

COMMON/VEX/XX ( 1 2 9 ) , Y Y ( 6 5 ) 

COMMON TIME , NI J 

COMMON NX , NY , NXl , NYl , PI , MAIN 

COMMON RATO , RATW , DELT , VELVAR , BXM , BYM , BZM , EYO 
COMMON I SEED , DY , DX , NIED , NI EDI , NSUMO 
C 

TP I = 2 . *PI 
DO 100 I = 1,2 
NSUMOl = NSUMO - 1 

I ONEX = I NT ( RAN ( I SEED ) * NSUMOl) + 1 

I ONPX = ION ( I ONEX , 1 ) 

IONPY = ION ( IONEX , 2 ) 

DO 10 J= IONEX , NI EDI 
DO 10 K = 1,5 

I ON ( J , K ) = ION ( J+l , K ) 

10 CONTINUE 

ION( NIED, 1 ) = I ONPX 
ION( NIED , 2 ) = IONPY 

RND = RAN ( I SEED ) 

DO 50 K=1 ,241 
DIFFl = FIW(K) - RND 
I F ( DI FFl .GE. 0.0) GO TO 55 
50 CONTINUE 

55 ION ( NIED , 3 ) = VXIW(K) 

VPAR = 0.0589* SQRT ( -LOG ( 1-RAN ( I SEED ) ) ) 

TPAR = TPI*RAN( ISEED) 

ION ( NIED , 4 ) = VPAR*SIN ( TPAR ) 

ION ( NIED , 5 ) = VPAR*COS ( TPAR ) 

NSUMO = NSUMO - 1 

100 CONTINUE 

RETURN 
END 

c 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE FOURN ( DATA , NN , NDIM , ISIGN) 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
REAL*8 WR , WI , WPR , WPI , WTEMP , THETA 
DIMENSION NN(NDIM) ,DATA( 16384 ) 

NTOT= 1 

DO 11 IDIM=1 , NDIM 
NTOT=NTOT*NN ( IDIM) 

11 CONTINUE 
NPREV=1 

DO 18 IDIM=1 , NDIM 
N=NN ( IDIM) 

NREM=NTOT/ ( N*NPREV ) 

IPl=2*NPREV 
IP2=IPl*N 
IP3=IP2*NREM 
I 2REV=1 



00 14 1 2 = 1 , I P2 , I Pi 

I F ( I 2 . LT . 1 2 REV ) THEN 

DO 13 11 = 12 , I2 + IP1-2 , 2 
DO 12 13=11 , IP3 , IP2 
I 3REV=I 2REV+I 3 — 1 2 
TEMPR=DATA( 13 ) 

TEMPI=DATA( 13 + 1 ) 

DATA( 1 3 ) =DATA ( I 3 REV ) 

DATA( 13+1 ) =DATA( I3REV+1 ) 

DATA( I 3REV ) =TEMPR 
DATA( I3REV+1 )=TEMPI 

12 CONTINUE 

13 CONTINUE 
END IF 

IBIT=IP2/2 

1 IF ( ( IBIT.GE.IPl) .AND. ( I2REV.GT.IBIT) ) THEN 

I 2REV= I 2REV-IBIT 
IBIT=IBIT/2 
GO TO 1 
END IF 

I 2REV= I 2REV+IBIT 

14 CONTINUE 
IFPl=IPl 

2 I F ( I FPl . LT . IP2 ) THEN 

I FP2 = 2 * I FPl 

THETA=ISIGN*6 . 2831853/( IFP2/IP1 ) 

WPR=-2 . *SIN( 0 . 5 *THETA ) * * 2 
WPI=SIN( THETA) 

WR=1 . 

WI = 0 . 

DO 17 13=1 , IFP1 , IP1 

DO 16 11=13 , I3+IP1-2 , 2 
DO 15 12=11 , IP3 , IFP2 
Kl = 1 2 

K2=Kl+I FPl 

TEMPR= ( WR ) * DATA ( K2 )-(WI ) *DATA( K2+1 ) 

TEMPI = ( WR ) * DATA ( K2 + 1 ) + (WI ) * DATA ( K2 ) 

DATA( K2 ) =DATA( Kl ) -TEMPR 

DATA( K2+1 ) =DATA ( Kl+1 ) -TEMPI 

DATA( Kl ) =DATA ( Kl ) +TEMPR 

DATA( Kl+1 ) =DATA ( Kl+1 )+TEMPI 

15 CONTINUE 

16 CONTINUE 
WTEMP=WR 

WR=WR*WPR— WI *WP I +WR 
WI =WI *WPR+WTEMP *WP I +WI 

17 CONTINUE 
IFPl=I FP2 

GO TO 2 
ENDIF 

NPREV=N*NPREV 

18 CONTINUE 
RETURN 
END 



